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SlMdARY 


A  two-dimensional  Lagrangian  finite-difference  computer  program  is 
developed  for  wave  propagation  analysis  of  buried  model  structures  under 
ground-shock  loading.  The  numerical  scheme  is  the  standard  method  originally 
proposed  by  von  Neuman  and  Richtmyer,  using  artificial  viscosity  to  smooth 
shock  fronts.  The  program  is  entirely  core-contained,  and  is  limited  to 
about  3,000  nodes  because  of  its  anticipated  application  on  personal 
computers  such  as  IBM-AT. 

Materia]  models  include  standard  hydrodynamic-elastic-plastic  relations 
as  well  as  a  new  equation  for  soils  and  concrete. 

Three  model  systems  were  considered  for  wave  analysis:  plane  slabs  with 
and  without  a  protective  soil  cover  and  a  buried  model  frame.  The  first  two 
represent  two  of  the  idealized  model  tests  described  in  Volume  I.  Since  few 
dynamic  data  exists  regarding  the  behaviors  of  the  sand  and  microconcrete 
used  in  the  construction  of  the  model  systems,  the  calculations  were 
intended  for  generating  the  qualitative  features  of  model  behaviors. 
Nevertheless,  the  computational  results  were  consistent  with  experimental 
observations  and  provided  a  rational  basis  for  interpreting  modes  of 
failure,  load  profiles  at  concrete  surface,  and  their  interrelationships. 

The  wave  analysis  of  the  buried  frame  indicated  that  modes  of 
structural  failures  under  dynamic  loading  can  be  predicted  by  directly 
focusing  on  shock  waves  that  excite  the  model  structure. 
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SECTION  I 


INTRODUCTION 


A.  OBJECTIVE 

This  phase  of  the  research  program,  undertaken  under  Contract  FO  8635 
85K  0052,  was  concerned  with  the  development  of  computational  capabilities, 
primarily  in  support  of  the  experimental  phase  of  the  program  described  in 
Volume  I  of  the  two-volume  final  report.  There  were  two  interrelated  goals. 
The  first  was  the  development  of  a  suitable  computer  program  and 
constitutive  equations  to  study  numerically  the  idealized  scaled  model 
testing  of  buried  structures  through  use  of  our  gun  facilities.  The  second 
was  the  demonstration  of  the  applicability  of  the  program  to  more  realistic 
model  structures.  Special  emphasis  was  placed  on  a  unified  treatment  of 
topics  that  are  traditionally  analyzed  separately.  They  are  ground «-  shock 
propagation,  soil-structure  interactions,  in-structure  shock  wave  propaga¬ 
tion,  and  failures  in  structural  members. 

B.  BACKGROUND 

Dynamic  problems  of  protective  structures  attendant  on  impact  and/or 
explosion  are  difficult  to  solve  because  of  nonlinear  differential  equations 
involving  complex  geometry  and  materials  behaviors.  Therefore,  traditional 
procedures  for  the  design  and  analysis  of  protective  structures  against 
conventional,  as  well  as  nuclear  weapons,  have  been  separated  into  two  or 
three  distinct  stages  (References  1  and  2).  For  instance,  in  the  case  of 
blast  or  ground-'shock  loading,  the  procedure  involves  two  separate  steps. 
The  first  is  the  determination  of  equivalent  loads  from  blast  or  ground 
shock  at  a  prescribed  standoff  distance.  The  second  is  the  design  analysis 
and  performance  predictions  of  a  given  structure,  based  upon  the  equivalent 
load.n,  by  assuming  an  equivalent  single-  or  multidegree  of  freedom  system. 
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Recent  improvements  in  the  development  of  predictive  and  analysis 
techniques  have  been  mostly  numerical  and  involve  computer  simulation  of  the 
dynamic  problems  by  use  of  finite  element  and  finite  difference  methods 
(References  3-6).  Large-scale  computations  have  come  into  widespread  use  in 
structural  and  ordnance  designs.  These  methods  have  now  been  developed  to  a 
point  where  they  can  handle  complex  shapes,  large  deformations,  and  failures 
involved  in  problems  of  impact  and  explosion.  The  principal  limitation  in 
the  use  of  these  numerical  methods  is  said  to  be  largely  the  uncertainty  in 
the  description  of  materials  behavior  (Reference  3).  If  the  response  model 
is  inadequate,  they  are  known  to  provide  solutions  which  are  not  even 
qualitatively  correct.  Therefore,  effective  use  of  any  of  the  numerical 
techniques  requires  iterative  adjustments  of  materia]  models  in  ('lose 
collaboration  with  material  and  structural  testings. 

The  present  investigation  leading  to  this  report  is  an  initial  attempt 
of  similar  type  to  analyze  and  aid  in  the  design  of  the  scaled  model  testing 
of  buried  concrete  structures  by  use  of  a  finite  difference  computet 
program. 

C .  SCOPE/APPROACH 

Many  computer  codes  are  available  fot  the  analysis  of  transient 
problems  associated  with  impact  and  explosions  (References  4-11).  Our 
compute)  program  is  a  two-dimensional  explicit  finite  difference  code  in 
plane  geometry  and  is  derived  from  other  codes  of  this  type  which  are  known 
by  their  acronyms  as  HEMP,  TWOODY,  STEALTH,  and  TROTT.  Special  features  of 
the  code  are  : 

•  the  code  is  core-contained, 

•the  cell  layout  is  easy, 

•  the  code  can  be  easily  modified  fot  new  materials 
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However,  it  is  restricted  in  the  size  of  problems  that  can  be  treated  and 
does  not  contain  features  such  as  slide  lines,  rezoning,  and  buffering  of 
cell  variables  that  are  common  in  large  general-purpose  codes  described 
above. 

Constitutive  models  in  the  program  are  presently  limited  to;  (1)  a 
standard  hydrodynamic-elastic-plastic  model  for  solids  and  (2)  a  model 
specially  developed  for  soils  and  concrete.  The  latter  is  developed  to  deal 
with  the  influence  of  porosity  on  the  inelastic  behavior  of  materials  in  a 
physically  consistent  manner  so  that  phenomena  such  as  shear  enhanced  pore 
compaction  can  be  represented.  New  models  will  be  added  as  needs  arise.  For 
instance,  a  high-explosive  equation  of  state  will  be  included  in  the  near 
future  to  consider  the  scaled  model  simulation  of  buried  structures  under 
close-in  detonations. 

Two  model  systems  were  investigated  numerically:  plane  slabs  with  and 
without  a  protective  soil  cover  and  a  buried  model  frame.  The  explosive 
loading  on  these  structures  was  simulated  by  a  shock-wave  loading  described 
in  Volume  I  of  the  final  repun.  Selected  results  from  the  numerical 
simulation  of  the  plane  slabs  were  compared  with  data  obtained  in  the 
experimental  phase  of  the  project.  They  are  a  qualitative  description  of 
failure  modes,  ground -shock  propagation  in  the  soil  cover,  and  structural 
response  behaviors  of  the  model  structures.  However,  because  of  limited 
experimental  data,  no  systematic  attempt  was  made  to  fit  the  experimental 
data  numerically  by  optimizing  material  parameters  in  constitutive 
equations.  The  principal  goal  of  the  calculations  was  to  demonstrate  the 
capabilities  of  the  code. 

The  remainder  of  this  report  is  divided  into  three  sections.  Section 
II  presents  a  summary  of  basic  continuum  mechanics  equations  for  describing 
two-dimensional  stress  wave  propagation  through  solid  and  porous  materials 
such  as  soil.  This  section  also  describes  the  algorithms  used  to  solve  the 
basic  equations  by  using  an  explicit  finite  difference  method.  Section  III 
describes  the  numerical  investigation  of  ground-shock  propagation  in  simple 
layered  systems:  plant'  slabs  with  and  without  protective  soil  cover  and  a 
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buried  frame.  Selected  results  from  the  first  example  were  compared  with 
experimental  results  reported  in  Volume  I  of  the  final  report.  Section  IV 
presents  the  conclusions  reached  from  the  sample  simulations. 
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SECTION  II 


CODE  DESCRIPTION 

A.  BASIC  EQUATIONS 

The  mathematical  description  of  transient  dynamical  problems  attendant 
on  impact  and  detonations  consists  of  the  conservation  laws  of  physics, 
initial  and  boundary  conditions,  and  material  models.  The  finite  difference 
technique  is  an  approximate  method  of  solving  these  equations  using 
discretized  time  and  space  coordinates.  The  Lagrangian  method  is  based  upon 
a  coordinate  system  that  moves  with  the  material.  It  has  been  developed  to 
a  point  where  general-purpose  codes  are  now  capable  of  handling  large-scale 
simulations  involving  50,000  nodes.  Since  there  are  many  excellent 
expositions  of  the  Lagrangian  technique  for  stress  wave  propagation  in  one-, 
two-,  and  even  three  dimensions  (References  7-11),  what  follows  is  a  summary 
description  of  the  basic  equations  and  computational  procedures  involved  in 
our  code.  For  details,  readers  are  referred  to  the  References  described 
above.  But,  on  exception  will  be  made  in  the  description  of  the  material 
model  for  soils  and  concrete  because  of  its  importance  in  the  solution  of 
problems  of  interest. 

1.  Conservation  Equations 

a.  Conservation  of  Mass 

The  equation  expressing  the  conservation  of  mass  is 

+  v,y  r  A/ A  (1) 

where  u  and  v  are  the  components  of  velocity  in  rectilinear  Lagrangian 
coordinates  (x,y),  A/A  is  the  aerial  strain  rate,  the  dot  signifies  the 
partial  time  derivative,  u,x  =  du/dx  and  u,y  =  dv/dy.  The  volumetric  strain 
rate  is  related  to  the  density  change  as  follows: 


-  ( ap/dt ) /p  =  u,x  +  v,y  +  W,z 

=  A/A  in  plane  x  -  y  geometry  (2) 

where  1/p  =  V  =  the  specific  volume,  and  K  is  the  velocity  component  in  the 
z-axis  (zero  in  plane  geometry). 

The  local  strain  rate  is  defined  as  the  symmetric  part  of  the 
velocity  gradient. 


f  =  u, 
xx  x 


t  =  v, 

yy  y 


£  =  (v,  +  u,  )/2. 

xy  ’  x  y 


The  local  rotation  rate  is  defined  by 


u  =  (v,  -  u,  )/2. 

xy  x  y 

b.  Conservation  of  Momentum 


The  equations  expressing  the  conservation  of  momentum  are 


pu  =  o  +  a  -  q,  (3) 

xx, x  xy,y  x 


pv-o+o-q,  (4) 

xy.y  yy.y  y  ' 


where  p  is  the  density,  u  and  v  are  the  components  of  acceleration,  <rxx, 
aXy,  and  o-yy  are  the  stress  components,  and  q  the  artificial  viscosity. 
The  sign  convention  for  the  stress  components  is  positive  for  tension. 
Pressures  are  positive  for  compression. 
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Inclusion  of  the  artificial  viscosity  tern  in  the  momentum 
equations  is  now  a  standard  technique  of  han  ig  discontinuous  shock  'waves 
numerically,  by  rendering  the  solution  continuous  using  viscous  effects. 
However,  care  is  necessary  so  that  the  viscous  term  does  not  affect  the 
solution  anywhere  except  at  shock  fronts.  In  our  code  the  main  artificial 
viscosity  term  consists  of  linear  and  quadratic  bulk  viscosity  components 
given  by 


q 


b ^ 2  p>t  +  b^A( p, 2/ p 
0 


for  r,^  >  0 
for  p,t  S  0 


(5) 


where  b]  and  b2  are  constants  and  Cs  is  the  speed  of  sound  given  by 


Cs  =  (3r/ap)g/2 


(6) 


where  s  indicates  the  differentiation  at  constant  entropy. 

When  a  large  distortion  of  computational  cells  is 
encountered,  additional  viscosities  are  included  in  the  deviatoric  stress 
components.  They  are  devised  to  minimize  hour-glass  shaped  distortion  of 
quadrilateral  cell  elements.  This  distortion  is  a  consequence  of  the  fact 
that  the  quadrilateral  cell  used  in  plane  geometry  has  eight  degrees  of 
freedom,  but  that  only  six  are  accounted  for  in  the  basic  equations  that 
provide  resistance  to  these  motions.  In  some  codes  these  viscosities  are 
known  as  the  triangular  artificial  viscosity.  These  triangular  stresses  are 
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where  b3  is  a  constant.  These  stresses  are  added  to  the  stresses  obtained 
from  the  constitutive  equations  described  in  the  next  section. 

c.  Conservation  of  Energy 


The  following  expression  of  the  conservation  of  energy 
ignores  thermal  heat  conduction.  This  is  a  reasonable  approximat ion  for  the 
time  scale  involved  in  stress  wave  propagation. 


E 


(p  +  q)  V  +  V(sxx 


e  +  s  £  4  s 

x*  yy  xy  zz 


£  4 

zz 


2  s  £  ) 

xy  xy 


(10) 


where  E  is  the  specific  internal  energy,  p  the  pressure,  sxx,  sXy,  and  syy 
are  the  deviatoric  stress  components  defined  by 


s  . 
ij 


a .  + 

1 J 


P 


(11) 


where  i  and  j  stand  for  x  and  y. 


2.  Constitutive  Equations 

The  present  program  contains  two  material  models.  They  are  a 
standard  hydrodynamic-el as tic-plastic  model  for  solids  and  a  special  model 
for  porous  media  such  as  soil  and  concrete.  The  former  model  was  initially 
developed  for  the  description  of  metallic  materials,  but  it  has  been 
effectively  used  even  for  materials  such  as  concrete  in  a  certain  range  of 
impact  conditions  (Reference  13). 

In  stress  wave  calculations,  it  may  be  necessary  to  permit 
fracture  of  materials  (say,  spalling)  during  the  calculations.  An  effective 
representation  of  such  a  separation  may  be  provided  by  letting  the  stress  in 
the  cells  along  one  side  to  reduce  to  zero.  The  algorithm  is  described  in 
Section  IIB. 
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a.  Hydrodynamic-elastic-plastic  model 

In  this  standard  model  the  stress  tensor  cjj  is  resolved  into 
the  pressure  and  deviatoric  stress  tensors  as  follows  (same  as  Equation 
(11). 


°i j  "  sij  P- 


The  hydrostatic  pressure  p  is  described  by  a  Mie-Gruneisen  equation  given  by 


p  r  yp)  +  p  r  e 


(12) 


where  r  is  the  Gruneisen  ratio  and  f^(p)  is  written  in  terms  of  a  polynomial 
function: 

f  (p)  =  a  (n  -  1)  +  a  (n  -  l)2  +  a  ( q  -  l)3  ,  (13) 

11  2  3 

where  i?  =  p/p0,  and  pc  =  initial  density. 


The  deviatoric  stress  components  account  for  elastic 
behaviors  and  are  calculated  by  the  frame- indifferent  isotropic  elasticity 
equations  (Reference  12), 
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ij  xj  ik  kj  kj  jk 


(15) 


V 

where  Sjj  is  the  co-rotational  stress  rate,  p,  shear  modulus,  is  the 

e 

component  of  the  rotation  tensor,  and  ejj  are  the  elastic  deviatoric  strain 
components  defined  by 


e . 
xj 


£  . 
XJ 
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(16) 
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n  plane  geometry  Equations  (14)  and  (15)  reduce  to 


®xx  ~  exx  +  2  *>  sXy 
S  yy  —  2(1  eyy  2l)  S  x  y 

®xy  =  2 (i c x y  ~  ®(sxx  -  syy)  (17) 


here  «  =  «xy  =  local  rotation  about  the  2-axis. 

The  transition  from  an  elastic  to  a  plastic  state  is 
ietermined  by  the  Von  Mises  yield  function 


g (a  .)  ^  s.  .s.  .  -  2YJ/3  .  (18) 
ij  ij  1J 

ihere  Y  is  the  yield  stress  in  simple  tension.  The  material  is  elastic  if 
{  <  0,  and  plastic  if  g  =  0,  whereas  the  condition  g  >  0  can  never  be 
realized.  In  the  plastic  state  it  is  assumed  that  the  total  strain  rate  is 
he  sum  of  the  elastic  and  plastic  strain  rates  and  that  the  plastic  strain 
ate  is  determined  by  an  associated  flow  rule  such  that 


€  .  . 
IJ 
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(20) 


In  this  formulation  the  plastic  strain  rate  becomes  normal  to  the  yield 
surface  as  expressed  by  Equation  (18).  This  relationship  provides  an 
axpedient  algorism  for  evaluating  stress  increments  in  plastic  state. 

b.  Model  for  Porous  Materials 


Continuum  plasticity  has  long  been  used  for  modeling  the 
mechanical  behaviors  of  geological  materials  such  as  soils  and  rocks 
(References  14  and  15).  Recent  models  involve  a  complex  combination  of 
mutiplastic  potential  surfaces  and  a  nonassociated  flow  rule  (References  16- 


10 


IS 


18).  However,  as  the  complexity  of  these  model  increases,  so  does  the 
difficulty  of  determining  their  material  parameters.  It  is  not  uncommon  to 
find  a  model  with  more  than  two  dozen  adjustable  parameters.  With  this  many 
parameters,  their  determination  is  rarely  complete,  particularly  when 
various  paths  are  involved  in  loading,  as  well  as  in  unloading. 

The  model  in  our  code  was  originally  proposed  by  Swegle 
(Reference  19)  as  an  extension  of  the  hydrodynamic  P-a  model  to  include 
shear  strength  in  the  description  of  porous  materials  including  geological 
materials.  The  most  important  feature  of  this  model  is  its  simplicity.  Our 
investigations  (References  20-21)  showed  that  as  few  as  two  to  four  free 
parameters  are  sufficient  to  deal  with  materials  such  as  metal  and  ceramic 
powders  as  well  as  various  types  of  soil.  Other  noteworthy  features  are: 
(1)  the  description  of  overall  stress  in  terms  of  the  stress  in  solid 
components  and  porosity,  and  (2)  the  use  of  associated  hardening  flow  rule 
to  describe  coupling  between  volumetric  and  deviatoric  inelastic  behavior. 

The  following  is  a  summary  of  Swegle’s  formalism  in 
incremental  form. 


a.  Stress 


Effective  stress  components  are  determined  by  those  of 
the  solid  components  and  the  solid  volume  fraction  a  as  follows: 


a  .  r  -p6 .  .  +  s. 
1J  1J  1J 


(21) 


p  -  otp  =  of  ( V  ,  E ) 
s  s 


(22) 


ij 


as 


ijs 


where  a  -  Vs/V,  the  subscript  "s"  stands  for  the  solid. 


(23) 
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b.  Strain 


Strain  components  are  partitioned  in  terms  of  volume 


components ,  i .  e .  , 


dr..  -  (dfl  /d0)d  r  .  . 
ijs  s'  1J 


(241 


where 


de  =  dV/V  and  d0  =  dV  /V  .  (25) 

s  s  s 

It  can  be  shown  that,  if  the  above  partition  of  the  strain  components  is 
used,  then  irrespective  of  deformation  modes, 

de. .  =  de.  .[1  +  (da/a) (dV/V)" 1 ]  .  (26) 

ijs  ijL 


where 


de .  .  =  or.  .  -  (d0/3)  6.  .  .  (27) 

ij  ij  ij 

c.  Elastic  regime 

Elastic  response  is  formulated  by  use  of  the  P-a  model 
and  a  frame  indiffernet  isotropic  Hooke’s  law.  That  is, 

(da/dP )eiastic  =  ^/Ko  ’  (28) 


where 


h(«)  =  1  ■»  [(1  -  a)  •  ( 1  -  «0)J(C0/Cso-  1)  (29) 

Kso  =  solid  bulk  modulus  at  zero  pressure, 

C0,CS0  -  sound  velocities  at  zero  pressure, 

dRijs  r  2Gsdeijs  +  ( ®iksk j~“k js  ik )^t >  (30) 
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d.  Plastic  Regime 

Plastic  state  of  the  material  is  determined  by  a  yield 
function  similar  to  that  of  the  perfectly  plastic  solid  described  in  the 
previous  section.  That  is, 

g  =  f(J},  J7  ,  a)  (31) 

where 

J1  r  °KK  "  _3p  ’  <32) 

and 

Jj  =  (1/2  s. .s.  .) 1/2  .  (33) 

N  2  U  U 

The  material  behaves  elastically  if 

g  5  0  (34) 

and  plastically  if 

g  =  0  .  (35) 

Then,  the  plastic  strain  increment  is  prescribed  by  the  associated  flow  rule 
such  that 


de  .  =  d .  +  dcp  .  (36) 

ij  ij  ij 

dcp.  =  d*  («g/*ff.  )  (37) 

where  the  superscript  describes  the  state  of  the  material. 


13 


There  are  many  potential  yield  functions  including  the 
well  known  Mohr-Coulomb  and  Drucker-Prager  models  for  the  description  of 
porous  materials.  Based  upon  our  experience  (References  20-21),  the  present 
code  contains  only  the  elliptic  yield  function  described  below.  But,  we 
also  tested  a  lemniscate  function  for  future  study.  Results  with  the  latter 
function  are  given  in  Appendix  A. 

e.  Elliptical  Yield  Function 

2  2 

g  =  [|p  -  ^  Pl^]  +  [(3  ^2)  /  Y](°)]  ~  1 

where 

Pm(«)  =  (Pl(«)  “  K( a) )  /  2 

P  (or)  =  (PP(«)  +  K(flt))  /  2 

K(a)  -  K0  +  cpp(o) 

Yi(a)  =  Qc( 1  -  a)1"  pp(o) 

Pp(«)  =  (y/fi)  [(1  ~  a)~2fi'/3  -  Cp  (1  -  a0)~2^^2\ 
and  c,  0o,  m,  Y,  0,  and  cp  are  constants- 

These  selections  have  been  made  based  upon  the 
observations  (References  20  and  21)  that  they  describe  the  general  features 
of  yield  surfaces  which  are  observed  experimentally  and  that  their 

parameters  can  be  understood  through  mechanistic  interpretation  of  flow 
mechanisms.  F01  example  the  current  form  of  pp(or),  which  is  the 

hydrostatic  inelastic  compression  of  porous  material,  is  determined  based 
upon  a  spherical  pore-collapse  model.  Detailed  descriptions  of  this  model 
are  found  in  Reference  20. 

3.  Boundary  and  Initial  Conditions 

Boundary  conditions  in  a  layered  system  may  be  divided  into  two 

distinct  categories:  external  and  internal  boundaries.  Typically,  external 


conditions  are  specified  in  terms  of  stress  or  velocity  components  for  all 
boundary  points.  Common  examples  are  a  free  surface  or  a  smooth  rigid  wall 
(zero  velocity  in  the  direction  normal  to  the  wall  surface).  Normally, 
boundary  conditions  at  an  internal  interface  are  the  standard  continuity 
conditions  for  normal  stress  and  velocity.  But,  they  could  become  very 
complex  if  interfacial  motions  such  as  sliding,  opening,  and  closing  are 
included.  In  the  current  code,  only  the  separation  of  an  interface  in  one 
spatial  direction  is  treated.  It  is  an  approximate  representation  through 
use  of  the  crack  opening  described  in  Section  IIB. 

Initial  conditions  must  be  given  for  all  dependent  variables  to 
determine  the  subsequent  motions.  That  is,  the  initial  state  of  stress  and 
strain,  as  well  as  geometry,  must  be  known. 

B -  COMPUTATIONAL  PROCEDURES 

In  a  finite  difference  approach  one  starts  with  the  governing 
differential  equations  and  approximates  them  by  appropriate  discrete 
equations  based  on  computational  grid  or  mesh.  Our  Lagrangian  explicit 
scheme  is  derived  basically  from  that  found  in  several  general-purpose  large 
codes  (References  8-12).  Therefore,  no  attempt  will  be  made  to  duplicate 
excellent  discussions  found  in  these  references.  We  will  limit  our  remarks 
to  a  description  of  the  inner  working  of  our  code. 

1.  Integration  of  Governing  Equations 

Figure  1  illustrates  our  finite-difference  grid  arranged  in  a 
staggered  rectangular  array  indexed  by  i  and  j.  Kinematic  quantities 
dealing  with  motion  such  as  positions  are  defined  at  integer  locations.  The 
remaining  quantities  such  as  stress,  strain,  and  internal  energy  are 
calculated  at  half-integer  points  as  averages  over  a  cell  volume  or  a 
surface  area.  For  example,  xjj  represents  a  Lagrangian  position  at  time, 
say  t  -  tn  and  Oj+1/2,  i+1/2  8  stress  component  for  the  quadrilateral  cell 
having  corners  at  (j,i),  (j+1),  i),  (j+1,  i+1),  and  (j,  i+1). 
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Computational 
cell  for 
a,  t ,  etc. 


Computational 
cell  for 
x,  x,  and  x 


M  j  j+1 

Lagrangian  Coordinate 

Figure  1.  Computational  Grid  Showing  Locations  of  Cells  for 

Computing  Kinematic  and  Equat ion-of-State  Variables 


The  integration  of  the  dynamic  equations  is  based  upon  the 
standard  leap-frog  method  and  proceeds  as  follows.  Initially  or  at  time 
t  --  tn,  all  quantities  are  defined  at  all  points  by  the  initial  conditions 
or  previous  calculations.  Computat ions  to  advance  one  time  step  from  tn 
to  tn^  is  done  in  four  steps,  as  illustrated  in  Figure  2.  First,  the 
momentum  equation,  Equations  3  and  4  are  solved  for  the  new  acceleration  at 
t  =  tn.  If  x*j  i  were  the  acceleration  being  calculated,  the  quantities 

required  are  specified  at  neighboring  four  half-integer  points.  These 
points  comprise  a  computational  cell  for  motion  variables  as  illustrated  in 
Figure  1.  Then,  the  resulting  acceleration  is  used  to  calculate  the  new 

n+i/?  n+1 

velocity  and  position,  xj^  xjj  by  using  time-centered  integrations. 


COMPUTATIONAL  SEQUENCE 


Figure  2.  Computational  Sequence  for  One  Time-Cycle 


1 


•  n+i/2  .n- 1/2  ,  ••  n  .,n 

x..  =  x  .  .  +  x .  .  At 

J.i  J.i  J.i 


n  -n+i/2  ..n+1/2 

=  x  .  .  +  x  .  .  At 
J.i  J.i 


where  Atn  -  tn+l  and  Atn+1/2  =  tn+1/2  -  tn  . 

The  order  in  which  computations  are  performed  is  such  that  at  the 
n+1 

time  the  position  of  Xj,i  is  computed,  the  positions  and  velocities  at  the 
other  vertices  of  the  quadrilateral  cell  having  smaller  j  and  i  are  already 
known.  A  schematic  of  the  order  is  shown  in  Figure  3. 


Points  already 
advanced  to  tn 
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Figure  3. 


0  Point  being 
advanced  to  tn+1 


o  Point  not  yet 
advanced  to  tn  +  1 


Order  of 
Computation 

Schematic  Representation  of  the  Order  of  Computation 
second  step,  we  use  those  new  positions  at  t  -  tn+l  to 


In  the  second  step,  we  use  those  new  positions  at  t  -  t n  J  to 
calculate  the  new  density  and  strain  components  from  the  equation  of  mass 
conservation  and  displacement-strain  relationships,  respectively.  These 
quantities  are  calculated  as  averages  over  a  computational  cell  shown  in 
Figure  1. 
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In  the  third  step,  using  the  new  density  and  strain  components, 
the  energy  equation  and  constitutive  equations  are  simultaneously  soivqd  for 
the  new  stress  components,  internal  energy,  and  other  constitutive  variables 
such  as  porosity.  Depending  upon  the  complexity  of  the  constitutive 
equations,  these  simultaneous  equations  may  involve  a  system  of  nonlinear 
coupled  ordinary  differential  equations  and  require  a  lengthy  iterative 
scheme  to  find  solution. 

The  fourth  step  is  a  preparation  for  the  next  time  cycle  and 
evaluates  the  next  time  increment  based  upon  the  Courant  stability 
condi t ion , 

At  -  tn"r]  -  tn  5  Ax/C  (41) 

s 

where  Ax  -  the  minimum  cell  dimension. 

In  our  program,  this  equat'r  is  modified  to  include  the  effect  of 
the  artificial  viscosity  as  fellows  (F-fereace  12). 

At  min  (Ax/C  )(1  -  31^)  (42) 

where  Cr  is  the  local  sound  speed. 

2.  Boundary  Condition: 

The  procedure  described  in  the  previous  section  applies  to  a  point 
in  the-  interior  of  the  homogeneous  Lagrangian  grid.  At  exterior  or  interior 
boundaries  th<-  algorisms  must  be  modified.  Currently,  only  a  limited  number 
of  boundary  conditions  is  provided  in  our  program.  lwo  kinds  of  exterior 
boundary  conditions  are  considered.  These  are  a  rigid  but  smooth  wall  and  a 
free  surface.  A  rigid  wall  (or  boundary)  is  represented  by  setting  the 
velocity  of  mesh  points  m  a  prescribed  direction  to  zero.  But  no  constraint 
is  imposed  on  the  motion  of  th‘j  mesh  points  in  the  direction  parallel  to  the 
wall.  At  a  stress  free  boundary,  since  the  stress  components  are  calculated 


as  averages  over  a  computational  cell,  a  massless  "phantom"  zone  is  created 
beyond  the  "real"  zone.  The  coordinates  of  the  outside  vertices  of  the 
phantom  cell  are  arbitrary  and  are  normally  set  equal  to  the  coordinates  of 
the  points  on  the  real  surface. 

For  interior  boundaries  the  present  version  of  the  code  cannot 
treat  sliding  internal  interfaces.  Therefore,  mesh  points  at  an  internal 
boundary  are  considered  to  be  common  to  both  sides  of  the  interface  and 
have  the  averaged  mass  of  the  adjacent  materials.  However,  an  effective 
representation  is  provided  for  a  separating  interface  in  one  direction  by 
letting  the  stress  in  the  cell  along  one  side  reduce  to  zero.  This 
procedure  can  also  be  used  to  represent  approximately  the  creation  of  a 
crack  in  the  interior  of  the  material.  The  behavior  of  cracked  material  is 
simulated  by  adjusting  the  stress  in  the  cell  so  that  there  is  no  norma] 
stress  across  the  crack.  But,  since  a  Lagrangian  cell  having  cracks  is  not 
allowed  to  separate  into  several  pieces,  the  stress  in  it  is  adjusted  to  the 
value  appropriate  to  such  cracked  material.  This  is  achieved  by  the 
following  stress  adjustments. 

o  -*  0 
xx 


o  ■*  a  +  AAr 

yy  yy 


o  ■*  a  +  Khi 

zz  Z7. 


a  -  0  (43) 

xy 


where  A  is  a  Lame’s  constant  and  hi  -  -  oxx/(K  +  2 jj)  . 


C.  INPUT  AND  OUTPUT 


Since  our  computer  program  is  a  special  purpose  code,  input  and  output 
data  are  kept  to  a  minimum  amount  necessary  for  solving  two-dimensional 
stress  wave  propagation  through  layered  systems.  The  input  data  consist  of 
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general  running  and  printing  instructions,  materials  data,  grid  layout  data, 
and  initial  conditions.  The  instructions  for  the  input  data  preparation  are 
described  in  Appendix  B. 


Because  of  our  emphasis  on  load  profiles,  the  current  output  are 
limited  to  stress  components,  positions,  and  particle  velocities.  However, 
the  program  can  be  easily  modified  to  print  other  variables  including  time 
histories  of  dependent  variables. 


SECTION  III 


SAMPLE  CALCULATIONS 

Two  types  of  problems  were  considered  for  sample  calculations  to  test 
the  computer  program.  The  first  is  the  calculation  of  stress  wave 
propagation  in  two  of  the  model  tests  described  in  Volume  T  of  the  final 
report.  ^he  second  is  the  analysis  of  the  model  testing  of  a  more 
realistic,  but  hypothetical  buried  structure.  However,  since  few  dynamic: 
data  exist  regarding  the  behavior  of  the  model  sand  and  microconcrete,  the 
calculations  were  intended  only  for  generating  the  qualitative  features  of 
model  behaviors. 

A.  SCALED  MODEL  SLABS 

Two  experiments  were  selected  for  sample  calculations.  They  are  M014 
and  M022.  The  major  goals  of  these  shots  were:  (1)  to  demonstrate  the 
capabilities  of  producing  a  selected  failure  mode,  say,  spalling,  by 
tailoring  shock  pulses  and  (2)  to  study  the  influence  of  a  protective  layer 
on  failure  modes  as  well  as  on  shock  profiles  (magnitude,  time  history, 
etc).  In  these  experiments,  the  dynamic  loading  was  modeled  by  impact  of  a 
projectile  fired  from  a  gas  gun  and  tailored  by  the  projectile  size,  speed, 
material,  mass,  etc.  An  idealized  test  configuration  used  for  the 
simulation  is  shown  in  Figure  4. 

Materials  properties  used  for  simulation  of  these  tests  are  summarized 
in  Table  1.  However,  since  no  systematic  experiments  were  conducted  for  the 
purpose  of  generating  dynamic  material  data  in  this  phase  of  the  research 
program,  no  attempt  was  made  to  optimize  these  properties  to  fit 
experimental  results.  Elastic  properties  are  those  found  in  standard 
handbooks.  Inelastic  properties  of  Lexan®,  aluminum,  and  microconcrete  were 
determined  using  the  von  Mise  criterion  and  stress  for  tensile  failure 
(Reference  22).  The  use  of  the  perfectly  plastic  model  for  the  concrete  was 
solely  due  to  the  lack  of  dynamic  data.  When  data  becomes  available,  it 
will  be  replaced  by  the  porous  model.  The  material  constants  for  the  sand 


TABLE  1.  MATERIAL  PROPERTIES  OF  ALUMINIUM,  LEXAN  AND  CONCRETE 


Aluminium 

Lexan 

Concrete1 

Bulk  Modulus  (Gpa) 

80 

3.47 

13.1 

Shear  Modulus  (Gpa) 

30 

0.90 

9.4 

Density  (kg/m3) 

2780 

1190 

1080 

Sound  Speed  (m/sec) 

6560 

1980 

3510 

Yield  Strength  (Mpa) 

75 

7 

28 

^Cunningham,  C.  H.  ,  Townsend,  F.  C.  ,  and  Fagundo,  F.  E.  , 

"The  Development  of  Micro-Concrete  for  Buried  Structrures , " 
University  of  Florida,  Gainesville,  198G. 


were  generated  through  use  of  Hugoniot  data  (Reference  23)  and  out  data  on 
wave  arrival  times  in  the  pressure  range  of  1-5  kb.  Figure  5  compares  of 
the  Hugoniot  data  with  those  gonerated  by  computer  simulation  of  shock  wave 
propagation  in  the  model  sand.  For  the  particle  velocity  of  less  than  0.2 
mm/  gsec,  waves  were  too  dispersed  to  define  a  meaningful  shock  front.  The 
calculated  Hugoniot  at  2  kb  was  fitted  to  that  estimated  from  the  wave 
arrival  times. 

1.  Shot  MO  14 

Figures  6  10  illustrate  selected  results  from  the  simulation  of 
Shot  MO  14  where  a  concrete  slab  was  shock  loaded  without  a  protective  sand 
layer  (see  Figure  4(a)).  The  aluminum  layer  was  used  to  prevent  impact 
damages  on  the  concrete  surface.  The  length  of  the  projectile  was 
arbitrarily  reduced  to  one  half  of  the  original  length  to  save  computer 
time.  This  change  has  no  influence  on  the  early-time  solutions  shown  in 
Figures  6-10.  In  this  simulation  the  separation  of  the  interface  between 
aluminum  and  concrete  was  provided  by  using  the  procedure  discussed  in 
Section  TIB.  Also,  the  edges  of  the  plates  were  assumed  to  be  free. 
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1.5  g/cc  Dry  Overton  Sand  —  1.5  g/cc  Dry  Overton  Sand 


Figure  5.  Comparison  of  the  Calculated  Hugoniot  Data  of  The  Model 
Sand  with  Those  of  Dry  Overton  Sand  (Reference  23) . 
Calculated  Shock-Wave  Velocities  Were  Fitted  to  Match 
Those  Measured  in  the  Model  Sand. 
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Figure  10.  History  of  at  the  Midpoint  of  Concrete 
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Figures  6  and  7  show  the  contour  maps  of  normal  stress  in  the 
longitudinal  and  transverse  directions  respectively.  The  progressive 
propagation  of  shock  wave  from  the  impacting  interface  is  typical  of  wave 
propagation  in  two  dimensions.  However,  the  stress  patterns  in  the  thin 
aluminum  plate  rapidly  become  very  complex  because  of  wave  reverberations 
between  the  two  interfaces:  Lexan®/aluminum  and  aluminum/concrete. 
Nevertheless,  as  shown  in  Figures  8  and  9,  the  stress  history  at  each 
interface  forms  an  expected  simple  "triangular"  profile.  Characteristic 
quantities  of  the  loading  profiles,  e.g.,  peak  stress  and  decay  time  are 
governed  not  only  by  two-dimensional  wave  interactions,  but  also  such 
parameters  as  impact  velocity,  impact  geometry,  elastic  constants,  etc.  The 
very  early  time  oscillation  in  the  history  of  normal  stress  at  the  second 
interface  is  a  normal  numerical  artifact  caused  by  a  sudden  increase  in  the 
rigidity  of  material,  i.e.,  concrete. 

The  appearance  of  strong  tension  in  the  concrete  slab  at  late 
times  in  Figures  6  and  10  is  a  wel 1-understood  phenomenon  and  was  caused  by 
the  reflection  of  a  triangular  compression  pulse  from  its  back  free  surface. 
When  the  tensile  stress  exceeds  the  ultimate  dynamic  strength,  a  fracture 
occurs  at  that  point.  If  the  fracture  extends  over  a  wide  region  as  shown 
in  Figure  6,  then  a  layer  of  material  may  even  split  away  (spall)  from  the 
rest  of  the  material.  In  localized  loading,  the  spalled  material  often 
takes  the  shape  of  a  cone  because  of  the  curvature  in  the  wave  front. 

Unfortunately,  since  our  fracture  model  does  not  describe  the 
process  of  dynamic  fracture  (Reference  24),  we  cannot  make  a  quantitative 
comparison  of  the  calculations  with  the  experimental  results.  But,  there  is 
a  qualitative  correlation  between  the  triangular  region  of  large  tension 
found  and  the  spalling  observed  in  the  test  (see  Figure  10  in  Volume  1). 

2.  Shot  M022 

Calculations  of  Shot  M022  are  shown  in  Figures  10-13.  The  major 
goal  of  this  test  was  to  examine  the  influence  of  a  sand  layer  on  loading 
profiles  and  the  modes  of  failure.  However,  in  the  simulation,  the  11 
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thickness  of  the  sand  layer  was  reduced  from  3-inches  to  1.18-inches  (3  cm) 
to  see  wave  propagation  in  the  concrete  plate  within  two  hours  of  computing 
time  on  our  IBM  370-168.  This  reduction  was  dictated  soley  by  limited 
computing  funds. 

Characteristic  features  of  wave  propagation  through  the  model  sand 
layer  are  illustrated  in  Figure  11  in  terms  of  the  histories  of  the  normal 
stress  cxx  at  several  successive  points  on  the  x-axis.  Two  noteworthy 
features  of  these  histories  are  the  rapid  attenuation  of  peak  stress  and  the 
dispersion  of  loading  profile.  In  general,  these  features  were  generated  by 
complex  two-dimensional  wave  propagation  through  the  model  system.  But,  the 
most  critical  parameter  is  the  slow  wave  speed  in  the  model  sand.  In  the 
pressure  range  of  our  interest  (less  than  5  kb)  this  speed  is  only  about 
one-  tenth  of  those  in  aluminum.  This  results  in  a  strong  lateral  unloading 
of  the  foward-moving  shock  from  the.  free  surfaces  of  the  projectile  and  the 
aluminum  plate. 

These  histories  indicate  that  by  the  time  a  wave  reaches  the 
surface  of  the  concrete  slab,  the  load  is  no  longer  a  triangular  shock 
pulse,  but  rather  a  step  load  having  a  time  scale  comparable  to  the 
fundamental  period  of  the  slab,  or  longer.  The  fundamental  period 
calculated  by  use  of  the  SAP  IV  program  (Reference  25)  was  about  140  fisec. 
This  means  that  the  sand  layer  very  effectively  transforms  a  highly 
localized  shock  pulse  into  a  long-time  structural  loading.  Therefore,  if 
there  were  any  failure  in  the  slab  it  will  be  one  of  the  structural  modes 
such  as  bending.  Figures  12  and  13  show  another  view  of  the  above 
described  transformation  in  terms  of  stress  contours. 

Although  the  wave  calculation  was  not  carried  out  long  enough  to 
observe  a  tensile  failure  caused  by  bending,  a  bending  failure  consistent 
with  the  calculation  was  observed  in  the  specimen  recovered  after  the  test. 

One  additional  conclusion  that  can  be  drawn  from  the  above 
described  results  is  that  the  wave  analysis  of  the  scale  model  testing  is  an 
effective  means  of  investigating  both  early-time  stress  wsve  response  end 
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"late-time"  structural  response  in  a  unified  fashion  by  directly  focusing  on 
shock  waves  that  excite  the  layered  system.  We  will  illustrate  this 
advantage  in  the  next  example. 

B.  A  BIT1IED  MODEL  FRAME 

The  goal  of  this  exercise  was  to  demonstrate  the  potential  use  of  the 
program  to  analyze  a  more  complex  model  than  a  plane  slab.  Of  particular 
interest  was  the  applicability  of  the  program  for  investigating  intermediate 
or  "long-time"  structural  response  in  a  model  structure  subjected  to  a  shock 
loading. 

The  geometry  of  the  problem  we  considered  is  shown  in  Figure  14.  It  is 
essentially  Shot  M022,  but  modified  so  that  the  model  structure  is  situated 
in  a  more  realistic  environment.  However,  the  selection  of  dimensions  was 
influenced  by  external  factors  such  as  the  cost  of  computing  time  and  the 
size  of  the  current  program.  The  soil  layer  was  further  reduced  to  1  cm. 

Since  early  time  responses'-  are  similar  to  those  found  for  Shot  M022, 
only  late  time  stress  contours  were  analyzed  to  focus  on  structural  response 
result  from  a  stress  pulse  propagated  through  a  layer  of  protective  soil. 

Figures  15  and  16  show  the  distributions  of  normal  stress  at  three 
successive  stages  of  wave  propagation.  Complex  profiles  in  the  aluminum 
plate  and  simple  wave  profiles  in  the  soil  layer  we'-e  very  similar  to  those 
found  for  Shot  M022.  In  spite  of  a  relatively  thii  oil  layer,  there  was  a 
subtantial  attenuation  and  dispersion  of  waves  c  id  mostly  by  lateral 
unloading.  As  in  the  case  of  M022,  the  unloading  w  determined  by  the  wave 
speed  in  soil  which  is  an  order  of  magnitude  slowe.  than  those  in  aluminum 
and  I.  ex  an®.  When  the  shock  pulse  reaches  the  concr  e  structure,  it  becomes 
a  slowly  varying  load  having  time  constants  comparable  to,  if  not  longer 
than,  1  he  period  of  the  first  structural  response  mode.  This  load,  in  the 
first  ordor  approximation,  is  a  distributed  normal  thrust.. 
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Therefore,  if  the  frame  is  excited  in  a  simple  structural  mode,  the 
mode  can  be  identified  by  examining  the  pattern  of  stress  distribution.  For 
instance,  it  is  believed  that  the  patterns  of  tension  and  compression  at 
45.7  |isec  (last  contour  maps)  correspond  to  the  shape  of  the  fourth 
fundamental  mode  (symmetric  bending  mode)  of  the  frame  shown.  The 
alternating  pattern  of  tension  and  compression  in  Figure  16  and  the 
corresponding  compression  regions  in  Figure  15  are  exactly  what  is  expected 
from  such  a  bending  response.  The  period  of  the  fourth  mode,  according  to 
the  SAP  IV  program,  was  about  53  jisec.  However,  no  attempt  was  made  to 
obtain  a  quantitative  comparison  because  of  the  difficulty  in  duplicating 
the  identical  dynamic  load  for  the  SAP  program. 

Nevertheless,  an  important  conclusion  of  this  calculation  is  that  model 
testing  combined  with  numerical  simulation  is  a  cost  effective  means  of 
generating  data  base  for  analysis  and  aid  in  the  design  of  protective 
structures  under  close-in  detonations. 
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SECTION  IV 


CONCLUSIONS 

In  the  foregoing,  we  described  (1)  the  successful  development  of  a  two- 
dimensional  Lagrangian  finite  difference  code  for  the  wave  analysis  of  the 
scale  model  testing  of  buried  structures  under  dynamic  loading  and  (2)  the 
development  of  constitutive  relations  for  the  description  of  porous 
materials  such  as  soil  and  concrete.  We  considered  three  sample 
calculations  involving  concrete  slabs  representing  actual  test 
configurations  and  a  buried  model  frame.  Results  of  the  former  agreed  with 
the  test  results  regarding  failure  modes,  and  provided  a  rational  basis  for 
interpreting  them  in  terms  of  load  profiles  at  concrete  structures. 

Some  of  the  advantages  of  the  wave  approach  are: 

•  regimes  of  both  stress-wave  response  and  early  structural  response 
can  be  analyzed  in  a  unified  scheme 

•  the  analysis  directly  focus  on  the  stress  waves  that  excite 
structures  as  well  as  shock-transmitting  layers 

Other  noteworthy  features  of  the  results  are: 

•  for  the  configurations  tested  here  the  slow  wavespeed  in  sand 
layers  is  the  major  cause  for  the  attenuation  and  dispersion  of 
ground  shock  that  changed  the  mode  of  failure  from  spalling  to  a 
structural  failure  by  bending. 

•  an  aluminum  plate,  because  of  its  fast  wavespeed  relative  to  sand, 
made  an  effective  layer  of  spreading  loads  over  a  wider  area  of 
sand.  This  phenomenon  may  be  of  some  interest  in  the  design  of 
protective  structures  by  using  the  concept  of  layered  systems. 
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APPENDIX  A 


GROUND  SHOCK  PROPAGATION 
IN  VARIOUS  SOIL  TYPES 

This  appendix  discusses  the  capability  of  the  porous  model  to  describe 
many  different  kinds  of  porous  materials  by  appropriate  choice  of  the 
functional  forms.  Presently,  the  model  involves  two  critical  functions  that 
determine  the  inelastic  behavior  of  porous  materials.  They  are  the  P-a 
relationship  and  the  yield  function  (see  Section  II. A).  However,  since  the 
present  P-a  relationship  is  determined  by  a  mechanistic  model  to  minimize 
the  number  of  free  parameters,  the  function  to  be  selected  is  only  the  yield 
f unct  i on . 


There  are  many  suggestions  for  the  yield  function  (some  models  involve 
more  than  one  yield  function)  to  describe  the  inelastic  behavior  of  porous 
materials  (Reference  15).  In  this  study  we  chose  a  lemniscate  function  and 
evaluated  the  predictive  capabilities  of  the  model  by  conducting  a 
parametric  simulation  of  the  scaled  data  on  ground  shock  propagation  in  five 
se:i  types  (Ref  er  ence  26'» .  Ttiese  soils  were  character  i zed  as:  (1)  loose 
density  sand,  (2)  medium  density  sand,  (3)  very  dense  sand,  (4.)  sandstone, 
and  (5)  silty  sand.  The  shock  data  is  said  to  be  a  compilation  of  more  than 
on*  hundered  explosion  tests  over-  tin  past  35  years. 

Tb*  behaviors  of  th<  five  soil  types  were  modeled  by  appropriate 
selections  of  the  material  parameters  in  the  lemniscate  function  illustrated 
in  Figure  A- 1  and  defined  by 

g  |(p  c())  /  A  ( a)  J  ‘  -*  [( /  B  ( a)  J  '  -  cos(r9  /  2 *)  (A.l) 


W*. 


A 


o 


}' 


a 


(  A  .  2  ; 
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DEVIATORIC  STRESS  INVARIANT 


LEMNISCATE  YEILD  FUNCTION 


9  = 


I  p*Cq 
\  A(cO 


PRESSURE 


Figure  A-l.  Lemniscatt-  Yield  Furction 


4K 


B(«)  =  c  cosB( »a  /  2)  •  A(a) 

D 

Pp(«)  =  (Y  /fi)  [d-a)_  2/,/  3  '  cp(l  -  aQ)"^/3] 

c0,  c^,  m,  /J,  and  Y  are  constants. 

The  specific  soil  properties  associated  with  the  yield  function  are 
summarized  in  Table  A-l.  The  relative  importance  of  different  parameters  is 
difficult  to  determine,  but  four  of  them  are  found  to  govern  characteristic 

features  of  the  ground  shock  propagation  in  the  five  soil  types.  They  are 

aQ  (or  the  initial  porosity  that  determines  the  sound  speed  ratio  Ct0/CS0,  n 
and  Cm  that  control  the  maximum  shear  strength  as  illustrated  in  Figure  A-l, 
and  ♦  that  determines  the  initial  slope  of  the  yield  function.  In  Suh’s 

formulation  (Reference  27),  ♦  is  the  slope  of  the  Mohr-Coulomb  failure 

surface. 

Other  material  properties  that  are  associated  with  the  Mie-Griineisen 
equation  for  the  solid  components  of  these  soils  are  assumed  to  be  the  same 
for  all  the  soil  types  and  are  listed  in  Table  A-2. 

The  parametric  simulations  of  buried  explosions  were  conducted  by  using 
the  gas  expansion  model  developed  by  Chadwich  et  al.  in  Reference  28  .  In 
this  model  the  pressure  of  the  exploding  cavity  is  given  by 


(A. 3) 

(A. 4) 


p  -  pQ(a/ao)  3r  (A. 5) 

where  pG  is  the  initial  pressure  when  the  cavity  radius  was  ac,  a  the 
current  radius,  and  r  a  constant.  The  explosive  parameters  were  obtained 
from  References  28  and  29  and  are  summarized  in  Table  A. 3.  Selected  results 
of  the  simulations  are  shown  in  Figures  A-2  -  A-5  and  Table  A-4.  In  these 
calculations,  however,  no  systematic  attempt  was  made  to  optimize  the 
materials  properties  to  obtain  the  best  fitting  to  the  experimental  results. 


49 


TABLE  A-l. 

PARAMETERS 

OF  THE 

LEMNISCATE 

YIELD  FUNCTION 

Parameters 

loose 

medium 

dense 

sand 

silty 

sand 

sand 

sand 

stone 

sand 

«o 

0.5 

0.59 

0.68 

0.78 

0.78 

n 

0.05 

0.05 

0.05 

0.05 

0.5 

♦(degree) 

53 

53 

53 

53 

45 

CG (MPa) 

0.  1 

0.1 

0.1 

0.  1 

0.5 

cm 

7.2 

6 

3 

1 

1 

Hi 

1.2 

1.2 

1.2 

1.2 

1.2 

Y(MPa) 

2 

2 

2 

2 

2 

A 

6 

G 

6 

6 

6 

CP 

0.95 

0.95 

0 . 95 

0.95 

0.95 

r.o 


# 

TABLE  A-2.  EQUATION  OF  STATE 


Propert ies 

loose 

medium 

dense 

sand 

silty 

sand 

sand 

sand 

stone 

sand 

Initial  porosity  (%) 

50 

41 

32 

22 

22 

Sound-speed  ratio  (Ct0/Cs0) 

0.05 

0.08 

0.15 

0.2 

0.2 

Density  of  Solid  (kg/m3) 

2700 

2700 

2700 

2700 

2700 

Shear  modulus  (Solid) (GPa) 

24 

24 

24 

24 

24 

Pulk  modulus  (Solid) (GPa) 

39 

39 

39 

39 

39 

*The  pressure-volume  equation  for 


sol  id : 


(V 

fV 

=  39 

S  -  1 

+  60.5 

-3-  -  1 

V 

V 

so 

so 

in  GPa 
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TABLE  A-3.  EXPLOSIVE  PARAMETERS 


Explosive  type: 

TNT 

» 

Charge  Weight,  L: 

512  kg 

w 

Radius  of  cavity,  au. 

A  2  cm 

Boundary  pressure,  PQ: 

22  Gpa 

• 

The  constant,  r: 

2.33 

Scaled  range  W: 

1  nr' kg 
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PEAK  STRESS  (MPa) 


Figure  A-2.  Attentuation  of  Peak  Stress  in  Various  Sands.  Lower  Curves 
Represent  Summary  of  the  Experimental  Data.  The  Bending 
of  Curve  D  at  About  50  MPa  Was  Caused  by  a  Too  High  Shear 
Strength  in  the  P-ct  Equation 


PEAK  VELOCITY  (m/sec) 


OCITY 


\  10.78  (L.S.) 

\  ; 

\\  0.78  (H.S  ) 


i  t  y  .  For  Clarity 
,  But  Kere  Located 
for  i  =  0.5,  0.59, 
\gain  the  Bending  of 
.ised  by  a  Too  High 


d/d 


PEAK  STRESS  DECAY 


P=  P, 


,e->t/ta 


W  =0.4075  (R- 3.26m) 


0  -J  1 - 1 - - - 1 - 1 _ - _ o  M-.o./  I 

0  1  2  3  4  5 


TIME  (x  10"3  sec) 


Figure  A-4.  Peak  Stress  Decay  in  lime  at  a  Fixed  Stand-Off 
Distance.  Experimental  Data  are  Only  Shown  for 
a  =  .78  (L.S.)  by  Broken  Lines 


TABLE  A-4.  CALCULATED  ATTENTUATION  COEFFICIENT,  n 


Sand  Type 

Tni t ial 

Peak  Stress  Level 

n 

«o 

(MPa) 

Loose  sand 

0.50 

1000  -  30 

3.52 

30  -  5 

2.80 

Medium  dense  sand 

0.59 

1000  -  60 

3.31 

60  -  7 

2.53 

Dense  sand 

0.68 

1000  -  60 

3.00 

60  -  10 

2.23 

Sand  stone 

0.78 

1000  -  100 

2.44 

(High  shear 

100  -  35 

2.05 

strength) 

35  -  20 

1.05 

Silty  sand 

0.78 

1000  -  80 

2.33 

(low  shear- 

80  -  20 

1.68 

strength) 

Figure  A-2  shows  a  comparision  of  calculated  attenuation  rates  of  peak 
stress  with  those  summarized  from  the  experimental  data.  The  calculated 
trend  of  the  attenuation  rates  is  in  agreement  with  that  discussed  in 
Reference  26,  involving  parameters  such  as  porosity  and  seismic  speed. 

A  similar  agreement  is  found  with  the  results  regarding  other 
parameters:  attenuation  coefficients  (Table  A-4),  attenuation  rates  of  peak 
particle  velocity  (Figure  A-2),  and  stress  decay  (Figure  A-3).  However,  no 
good  correlation  was  obtained  with  particle  velocity  decay  (Figure  A-4).  In 
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:hese  figures,  the  stress  and  velocity  decays  that  characterize  pulse 
•rofiles  are  defined  by  Drake  and  Little  as  follows  (Reference  26).  First, 
:he  arrival  time  ta  is  defined  by 

t  a  =  R/c  ( A . 6 ) 

where  R  is  the  distance  from  the  explosion  and  c  is  the  seismic  wave 
propagation  velocity.  Then,  the  stress  and  velocity  decays  are  defined  by 

0(0  =  po  exp(-rt/t  )  (A. 7) 

v ( t )  =  v  (1  -  /Jt/t  )  exp  (-/Jt./t  )  (A. 8) 

O  a  a 

where  ap  and  vp  are  the  values  of  the  peak  stress  and  peak  particle  velocity 
and  y  and  p  are  time  constants.  They  find  that  for  most  applications  these 
time  constants  may  be  approximated  by  r  =  1.0  and  p  -  1/2.5. 

They  also  found  that  the  rise  time  of  these  wave  forms  is  typically 
about  one/tenth  of  the  travel  time  to  the  target  point.  The  corresponding 
calculated  results  varied  from  0.06  to  0.17  depending  upon  soil  types  and 
locations.  An  i ccurate  evaluation  of  the  coefficient  is  difficult  because 
of  the  ambiguity  in  the  definition  of  travel  time  for  a  dispersed  wave 
profile.  But,  the  calculated  results  are  consistent  with  the  reported 
empirical  value. 
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APPENDIX  B 


INPUT  PREPARATION 

Four  sets  of  input  data  are  required  to  run  a  problem.  They  are,  in 
the  order  of  appearance,  general  running  data,  materials  data,  cell  and 
coordinate  layout  data,  and  velocity  data.  As  an  example  the  complete  input 
deck  for  the  calculation  of  the  buried  frame  is  shown  in  Figure  B-l.  Any 
consistent  system  of  units  can  be  used  for  the  calculations.  Some  popular 
systems  are  shown  in  Figure  B-2  (Reference  9). 

A.  RUNNING  DATA 

Line  1  (215,  F10.0 
Columns 
1  -  5 

6-10 
11-20 

Line  2  (515) 

Columns 

]  -  S 

-1  -  Fixed  y-velocity  at  i  =  imin  and  imax 

-2  =  Fixed  y-velocity  at  i  =  irain  only 

-3  =  All  free  edges 

-4  =  Fixed  y-velocity  at  i  =  imin> 

Fixed  x-velocity  at  j  =  and  jmax 

-5  =  Fixed  y-velocity  at  i  =  imin  and  i, 

Fixed  x-velocity  at  j  -  jmin  • 


Variable 
boundary  condition 


Variable 

frequency  in  cycles  for  printing  selected 
solutions 

the  maximum  number  of  computing  cycles 
initial  time  increment 


C  -  10 
11  -  15 


number  of  blocks  in  the  problem 
number  of  materials 


30 

450 

.  IE-06 

* 

General 

-2 

7 

4  1 

running 

. 3470E+ I } 

• 

dat  a 

119  E~ 

01 

. 689E08 

. 904E10 

.  198F0C. 

5 

4. 

.  1 

0.02 

300000. 

.800 

E12 

.278 

E01 

.  750E09 

. 300E12 

-  656200 

4. 

.  1 

0.02 

Mat  erials 

3 . 900E1 1 

6.05E1 1 

1 

i 

data 

2.7 

7.5  E08 

2.40E11 

512799. 

4. 

.  1 

0.02 

.08 

.59 

. 131E12 

. 208F01 

. 282E09 

. 940E1 1 

. 351E06 

4. 

.  1 

0.02 

1 

15 

0. 

2.8 

2.8 

0. 

1  | 

1 

11 

0.0 

0.0 

2.0 

2.0 

j 

15 

1C 

2.8 

3.0 

3.0 

2.8 

1 

1 

11 

0.0 

0.0 

2.0 

2.0 

j 

10 

21 

3.0 

3.7 

3.7 

3.0 

2 

1 

41 

0. 

0. 

8.0 

8.0 

Cel  1  and 

21 

31 

3.7 

4.7 

4.7 

3.  7 

3  | 

coordinate 

i 

X 

41 

0. 

0. 

8.0 

8.0 

1 ayout 

21 

61 

4.7 

7.7 

7.7 

4.7 

3 

I 

20 

41 

5.0 

5.0 

8.0 

8.0 

31 

61 

4.7 

7 . 7 

r>  r-> 

t  .  / 

4.7 

‘1  ; 

21 

2C 

4.0 

4.0 

5.0 

5.0 

31 

41 

4.7 

5.7 

5.7 

4.7 

4 

1 

2 1 

0.0 

0.0 

4.0 

4.0 

l  i 

if; 

15000. 

"1  Velofits  d 

Figure  B-  1 . 


Input  Data  for  the  Calculation  of  the  Buried  Franu. 
Shown  in  Figure  14 


f.O 


Sets  of  Units 


Quantity 

c.g.s. 

c.g.  ns 

S.  I. 

f.p.s. 

Time 

s 

ns 

s 

s 

Length 

cm 

cm 

m 

ft 

Mass 

gm 

gm 

kg 

slug 

Force 

dyn 

T  dyn 

Newton 

lb 

Energy 

erg 

T  erg 

Joule 

ft.  lb 

Energy  Density 

erg/gm 

Mbar  cm3/gm 

J/kg 

ft . lb/slug 

Power 

erg/s 

T  erg/s 

Watt 

ft. lb/s 

Densi ty 

gte/cm  3 

gm/cm  3 

kg/m  3 

slug/ft  3 

Pressure 

dyn/cm  2 

Mbar  Pa 

lb/ft  2 

Figure  B-2. 

Consistent  Systems  of  Units  That 
Code  Calculations 

Can  Be 

Used  for 
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16-20 


velocity  initialization 
1  =  velocity  initialized  for  all  i  up  to 
an  interface  j  value  for  a  projectile 
impact 

-1  =  velocity  initialized  for  all  i  and  j  from 
an  interface  value  to  jma>:. 

21  -  25  maximum  value  of  j  in  initial  calculations. 

P.  MATERIAL  DATA 


Line  1  (4F10.0,  15) 


Columns 

Variables 

1  -  10 

aj  in  Equation  (14) 

11-20 

&2  in  Equation  (14) 

21  30 

arj  in  Equation  (14) 

31  -  40 

r  in  Equation  (14) 

1 

Cl 

1  -  porous  material 

Line  2  (4F10.0,  15) 

Columns 

Variables 

1  -  10 

specific  density 

11  -  20 

yield  stress 

21  -  30 

shear  modulus 

31  40 

longitudinal  sound  speed 

41  45 

fracture  indicat  or 

5  r  fracture  in  the  x-direction 

6  -  fracture  in  the  y-direction 

Line  3  (3F10.0) 

Co) umns 

Variables 

1  10 

quadratic  viscosity  constant 

11-20 

linear  viscosity  constant 

21-30 

triangle  viscosity  constant 
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Line  4  (2F10.0: 
Columns 
1  -  10 

n  -  20 


Vari ables 

^to  /  CSo  =  sound  speed  of  porous  material  ' 

sound  speed  of  its  solid  component 


Line  5  (FIO.O; 

Columns 
1  -  10 


skip  this  line  if  the  fracture  is  not 
considered. 

Variables 
fracture  strength 


Repeat  Line  1-5  as  many  times  as  the  number  of  materials  in  the 
problem. 

C.  GRID  LAYOUT  DATA 

Line  1  (2I5.4F10 . 0 .  15) 


Columns 


Variables 


1 

-  5 

jl 

6 

10 

j2 

11 

-  20 

XA(  1 ) 

21 

-  30 

XA(2) 

31 

-  40 

XA(3) 

41 

-  50 

XA(4) 

51 

-  55 

material  number 

where  jl  and  j2  define  Lagrangian  x-posilions  of  a  quadrilateral  block 
as  shown  in  Figure  B.3.  XA(n)  are  x-coordinates  of  the  corners  of  the 
block  read  in  counterclockwise  direction  starting  with  pint  of  smallest 
j,  i  values. 
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Line  2  (215.  4F10.0) 


Columns  Variables 


1 

-  5 

il 

6 

-  10 

i2 

11 

-  20 

YA(  1) 

21 

-  30 

YA  ( 2 ) 

31 

-  40 

YA(3) 

41 

-  50 

YA(4) 

where  il  and  i2  are  Lagrangian  y-positions  defined  in  Figure  B-2  and 
YA(n)  are  y-eoordinates  of  the  corners  of  the  quadrilateral  block. 

C.  VELOCITY  DATA 


Line  1  (315.  F10.0) 


Columns 
1  -  5 
6  -  10 
11  -  15 
16-25 


Variables 

minimum  i  value  initialized  at  the  velocity  u 
maximum  i  value  initialized  at  the  velocity  u 
interface  j  value  for  a  projectile  impact 
initial  velocity  of  a  projectile 
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APPENDIX  C 


LISTING  OF  COMPUTER  PROGRAMS 


6G 


IBPLICIT  REAL‘3  (A-B.O-Z) 

CCB^CN/T/X  (3000)  ,X  (3000)  ,  ID  (3  000)  ,XD(3000)  ,8(3000)  ,  A  (3000) 

.  Z (jOOO) , D(3000) ,  5XX (3000) , SX  I (3000 ) . SZZ (3000) , TXT (3000) 

♦  , TXX  (30  00)  ,  TXX  (3000)  . TZZ  13000 ), P  (3000 )  ,£(3000)  ,  XX  13000)  , 

.LVV  (3000)  ,  EXX  (3000)  ,EXX(3000)  ,EZZ{3000)  ,EXX  (3000)  #  ALP  (3000)  , 

♦  Id  (JO 00) 

COS.ION/EWS/tuSTC  (to)  ,  EWSTD  (to)  ,  E'jSTS  (6)  ,BHO  (6)  ,  XC  (6)  ,  UBO  (6) 

♦  ,CLIMt),C.Su(6),TfcI0(b).SP(6),EwSTG(6)  „ 

COBBG  B/GER/UZ  ZSO,DT,DTW,DTW,IJ8tJRD,JBAX,JBIN,KI1AX,Kfl.IN,KCHEK 
CCM BC  R/ 1 B  0/  BPS  (10) 

CGaBOR/PSF>/  TSfi{10) 

COBBOB/POR/  ALP (3000)  ,SXIS  (30  00) ,  STXS (3 00  0)  ,SZZS(3000) ,  TXXS (3000) 

♦  ,  TT  X  S  (3000)  ,  TZZ  S  (3000)  ,TX  IS(3000)  ,?S(3000)  ,EVVS(3000)  ,EXXS{3000) 

», EX  IS  (3  000)  , EZZ  S ( 3000 ) , EXXS (3  000)  ,EJ(3000) , ALPO  (6) , I  POE (6) ,CTCS{6) 

DIBEBSION  XI  (  100,  100)  .XL  (100,  100)  ,HH<  100,  100)  ,LVAB{100,  100) 

J  SI ZE  =  1 5000 
JIX= 100 
RXX=  ICO 

RE AO  (  1,2222)  NB.NBAX.DT 
PCSBAT  (215, £10.  3) 

•  BITE  (3,2223)  KK.KBAX.DT 

**  Oj  FOP.  BAT  (  IX,  *KH='  ,13,31,'  KBAX='  ,15,  3X,'DT='  ,E  12.5) 

c 

CALL  LAYOUT  |J  SI  Z E,  J  aX  ,K  XX, XL,  XL ,  B H,  LV  AR) 

T  I  H  E  =  3. 

DT  S  •=  0  T 
U-  Q 

BAN- N  B 
IvO  N=«*1 

IF  (0. GT.UBAX)  CALL  EXIT 

CALL  SWEEP  (JSIZE.JXX.KXX. X~, XL, BB.LVAB) 

L»  1  =  LVAB  (  1  I  ,2) 

LV2=_VAB  ( \o,2> 

LV3^LVA3  (2  1,2) 

LVu  =  LVAB  (26,2) 

LV5  =  LVAii(31,  0) 

TIf.E=TISE*DI 

CC-  WRITE  (3,  103  1)  N.TIBE, P(LV1J  ,EV?(LV1)  ,ALP(LV1)  ,  SIX  (LVl) 

U-w  •  ,  S  j  X  (LV  I)  ,  P  (LV2  )  ,  E  V  V  (LV2)  ,ALP  (LV2)  ,  SXX  (LV2)  ,SXX  (L  V  2) 

cc-  B  BITE  (9  ,  1030)  N  ,  TI  BE.  TXX  (L  V  1)  ,TXX  (LV2)  ,  TX  X  (LV  3)  ,TXX  (LV4)  ,  TX  I (LV5 ) 

.  WrITE  (9, 103  0)  H  .  TIHE ,  TXX  ( LV  1)  ,TXX|LV2)  ,  TX  I  ( LV3)  ,TXX  (LV4)  ,TXX  (LV5) 

WRITE  (8,  1030)  N,TIf.E,P  (LV  1)  ,  E  V  V  (LV  1 )  ,  ALP  (LV  1)  ,  X  0  ( LV  1)  ,  E  J  ( L?  1) 

.RITE  (9  ,  1  0  3  0)  N.Tir.E,  TXX  (LV  1)  ,TXX(LV2)  ,TXX(LV3)  ,TXX  (LV4)  ,IXX(LV5) 
K  3  1  P06BAT (15,  1  IE  11.4) 

C - TAPi  o  FOB  M, II  BE, DT 

>_  . RITE  (8,8 1 1 )  B, T 1  HE , DT 

C  FCSBAI  (15, 2E  10.  3) 

IF  (k.  BE . BUM )  GO  TO  951 
W  SI TE  (3,2224) 

FOBBAT  (3 X , ' J ' , J X,  •K,,2X,,LVAEH', 1 IX, 'X*  ,  lOX.'XD', 9X, 'TXX* , 
*8X,'TXX',12X,,P,,9X,"TXT',9X, 'SXI» ,9X, ' SXX’, 8X, 'ALFA* .//) 
RBR-BBN+NR 

IF  (W. LT.  145)  GO  TO  95  1 
DC  950  K= 1 , KCHEK 
DO  920  J=  1, JfllX 


L  V  Afifl-L  V A  B  (A, J) 

IP  (LVARB.LE.O)  GO  TO  920 

IP(fl  (LVABB).LE.O. ABD.O. EO.  1)  GO  TO  920 

WBITE  (3,  1040)  J  ,  A,  LV  ABB  ,  I  (L  V  A  B  B)  ,  ID  ( L  V  A  BH)  ,  TX  X  (  L  V  A  B  H) 

♦  .TXX  (LV  AEH)  ,P  (LVABB) , TX X (LV A£ B)  ,SXX  (LV ABB) , SI X  (LV ABB)  ,ALP  (LVABB) 
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C - DISK  9  FOB  N,J,K,I,T,XD,TD,P,E 

c  DISK  10  POE  K,J,K,EXX,EYI,EIY,EZZ,ALFA,EVV 

C  DISK11  FOB  N,  J,  K.TXX.TY  Y.TIt.TZZ,  SIX,  STI.  SZ2 

v 

C  N  BIT  E  (9 , 9  1  I )  N,  J,  K, X(LVARfl)  ,T  (LVAEH) ,XD  (LVABH) , ID  (L7AHB) 

<-  ♦  , P  (LYABfi)  ,1  (LVABH) 

<-  *11  FOBBAT  (319, 6E  10.  3} 

C 

C  BELTS  (10,9  I  1)  N,  J.A.EXI  (LVABH)  ,  EYY  (LVABH)  .EXT  (LVABH)  , 

C  ♦  EZZ  (LVABH)  , ALP  A  (LVAEH) ,EVV  (LV AEK) 

B File  .  111  1)  h.J.K.l  (LVABH).  X  (LVABH)  ,TXX  (LVABB) , TY I  (LVABH)  , 

♦  TZZ  (LVAhBj  , 1X1  (LVAEH) 

1  I  11  FOEBAl  (31 5, 6E  10. 3) 

U.  HEITE  (10,  111  I)  N,J,K,X(LVARH)  ,  Y  (LVABB)  ,  XD  (LVAEH)  ,  P(LVABH)  , 
CLC  *TXX 'LVABH)  ,TYX (LVAEH) 

*20  CC'  '  I  N  ;J  E 
) SO  CONTINUE 
49  1  CONTI  HUE 


‘j  T  N  -  O  I 

DT- SHIN  1  (0.  9*DTB,  DflAX  1  (1.  2*DT ,0. 0  3S»DTV)  ) 
DTN  =  0 . 5* (DT  ♦  DTN) 

IF  (LT.  31.  1.  E-  12)  GO  TO  100 

K  HITE  (3  .1020) 

1020  FCBBAT  (  IX, • STABILITY* ) 

10  30  FOEHAT (15,7110.  3) 

13-0  r OHBAT (315, 21 , 10E  12. 5) 


ST  0  f 
EM 


S  U&iOUT I NE  LAYOUT  (JSI ZE , J XX .1 XX , X L, Y L , H H,  LVAB) 

1HPLICIT  BEAL*8 (A-H.O-Z) 

ConaoN/T/X (3000) , X (3000) ,XD  (3  000)  ,ID (30  00) ,H(3000) , A (3000) 
♦,Z(3000)«O(3000) , SIX (3000) ,SII{3000) ,SZZ{3000) ,TXY(3000) 

♦  «TXX  (3000) ,TYY (3000) , TZZ  (3000) , P { 3000) , E ( 3000)  , YT (3000)  , 

♦EVV  (3000)  ,  E  X  X  (3000)  ,EII{3000)  ,EZZ  (3000)  ,EXY  (3000)  ,  ALP  (3000)  , 

♦IH (3000) 

CCBflON/EkS/lu3TC(6)  ,  EwOTD  (6)  ,  EQSTS  (6)  ,HHO  (6)  ,  IC  (6)  ,  UBU  (6) 

♦  .CLIN  (o)  ,CVS0  (6)  ,TBIU  (6)  ,SP  (6)  ,£QSTG(6) 

CGHHON/GE  Si/UZEBO,  DT  ,  DTN  ,DTW,IJBUND,JHAX,JflIN,  KB  AX  .KBIN.KCUEK 
CCHBON/IND/  NFR(10) 

COHHON/FSE/  T SB  (10) 

COflflON/POB/  ALP (3000)  #SXXS  (3000) , SYYS (3000)  .SZZS (3000) ,TXIS  (3000) 

♦  ,TYIS  (3000)  ,TZZS (3000) ,TXYS  (3000) , PS  (3000) ,BVVS  (3000)  , Eli S  (3000) 

♦  , EYYS  (3000)  , EZZ  S (3000) , EXYS  (300  0)  ,EJ  (30O0) , ALPO(6) ,IPOB (6) ,CTCS (6) 
DIMENSION  XL (KXX , JXX)  ,IL(KXX, JXX)  ,BB(KXX,JXX)  , L V A B  (KX X , JX X) 

♦,XA(4),YA{4) 


L -  10.  J=  I,  3000 

x  (j  j  -  o. 
x  (J>  -■  0. 

ID  (0)  =0. 

YD  (J) -0. 
a  {Jj  =  0. 

A  (J)  =  0. 

.  (J  i  =C. 

D  ( J ;  -  0. 

SXX  (J)  =  .. 

SYY  ( J  )  =  0. 

SZ-  (J)  =0. 

I  XT  (0  )  -  0. 

TXX  ( J )  =  C . 

T1Y  (J)  -  0. 

TZZ  (0)  -  0. 
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P  (J)  =0. 

£  (J)  =  0. 
yy  (J) =3. 


IIP  (J)  =  0. 
SXXS  (J) =0. 
SYYS  (J) =0. 
SL2S  (J) =0. 
TXXS  (J) =  0- 
TYYS  (J) =0. 
T2ZS  ,J) =3. 
TIYS  (J)  =3. 
PS (J) =0. 

EV  VS  (J)  =0. 
EXXS ( J j =0. 
ETYS (J) =0. 
E2ZS  (J) =0. 
IXYS (J) =0. 
EJ  <J)  =3. 

10-  CCATIN'JE 


J  K  =  JXX»KXX 

Dc  104  1=1, EXX 
DC  10*4  J  =  1,  JXX 
XL (I. J) =- 9  y  S. 

XL  ( I  .  J) 

MH  (I  ,  J)  =  0 

10-  LVAP  (I, J) =0 

KHAX-KXX 
JHAX= JXX 

LEAD(I,IOIO)  IJbUHL. N BLOCK, MM Tfa S, IV TYPE, KCHEK 
MBIT E  (3,2  1.  1)  10  BUKO,  KBLOCK,HtiTRS,IVTYPE,  KCHEK 
2  121  FOB  HAT  (  lx, ■ IJ3UHD=' , 12, 3X , ' MBLOCK=* ,12, 3X,* NHTB  S  = • , X2.3X, 

♦ ' 1VTYPS=', 12, 3X , ' KC  8EK= ',!-,//) 

OG  10  5  i=  1,  M.-.7SS 

EE  AO  1  I,  10  11)  Ek*STC  (1)  ,E*STD  (I  )  ,  EQ5TG  ( I)  .EgSTS  {I )  ,  XPOP,  (I) 

M&ITt  (3,2  12  2;  E*STC(1)  ,  E  jSTO(I)  ,EQ3TG  (1)  ,  EgSTS  (I)  ,IPOB(I) 

- F0SflAT <  IX, •E*STC=', El  2. 5, 3X, ' EjSTD=*,E12.  5, 3X , • EgSTG= • ,  E 1 2.  5 

♦,3X,*EgSTS=\E12.5,3X,‘IPOR  =  ',I2,/) 

BiAO(l,  1312)  hnO(I)  ,YC(I)  ,3 HU  (I),SP(I)  ,  NF3(I) 

WSITE  (3,212  2)  BhO(I)  ,  YC  (I)  ,UMU  (1)  ,3P(I)  ,»F8  (I) 

-12  3  FOB  HAT ( IX, • EE NS  IT Y= • , El  2. 5, 3X , *10=' , E 12. 5,3X, 'SHE IE  HOD.=  •, E 1 2. 5, 
♦ 3X  ,  ' 5  OU  HD  SPEED =' ,t 1 2 . 5 , 3  X , ' F  B  A  CT  U  3  E  0P7.  =  »,I2,/) 

BEAD (  1,  1  104)  CJS3(I)  ,CLIH(Ii.TBIC(X) 

•  EiTE  (3,221.)  CgSg(I)  ,  C  LI  M  (I)  ,TBIg(I) 

,^1S  FOBHXT  (IX, • eg sa  =  ',r  10.5, 3X,'CLIM= '.FID. 5, 3X,'TBIQ=< ,P10.5,/) 
IFJIPOfi  (I). EQ.O)  GO  TO  114 
E  E  A D  (  1 ,  1011)  CTCS  (I)  ,  ALPO  (I) 

MB  IT  E  (3,323  3)  CTCS  (I)  ,  ALPO  (I) 

3^33  FOB  HAT  ( IX , • CTCS  = • , E 10 . 3 ,31 , • A LPO= • , F 10.  5,//) 

11-  IF  (HFB  (I)  .EC.  3)  GO  TO  105 
BEAD (  I,  10  1  1}  73 B  (I) 

■  SITE  (3,4344)  TSS  (I) 

-3  44  FORMAT <  IX, ' FEACTUBE  STB ESS= ' , E I  2. 5) 

105  CCNTIMUE 


C--CELL  LAYOUT 

V- 

DO  250  HB=  1 ,MBLOCR 

BEAD  (  1,  1030)  K 1,5 2,  (X A  (I) ,1  =  1,4) , HAT 
HBITE  (3,6580)  K  1 , K2,  (XA  (I) , 1=  1 , 4)  , H AT 
B  BAD  (  1,  1030 )  J1,J2,  (IA(I)  ,1=1,4) 

HBITE  (3,6590)  J  1,J2,  (YA  (I)  ,1=  1,4) 
e-SC  F0BHAT(1X,,K=',3X,2I5,3I,'X=',4E12. 5, 3X, • HAT= ' , 13) 
t,590  FOBHAT(IX,'J='.3X,2I5,3X,'r  =  ',4E12.5,/) 
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1) JDK=  (J2-J1)  •  (K2-K1) 


DO  210  K=  K 1  ,  K2 
DO  210  J=J1,J2 

If  (XL  (K,J)  .  £J. -999. )  XL|K,J)=  (  (XA  (1)  *  (J2-J)  *14(4)  •  (J-J1)  )  *  (K2 

♦  K)  ♦  (X  1  (2)  •  (  J2-J)  *Xk  (3)  *  (J-J1)  J  •  (K-K  1)  )/DJDK 

If  (XL  (K,J)  .  IQ.  -999.)  IL(K,J)=  (  (XA  (1) *  (J2-J)  *14(4)  *(J-J  1))  *(K2 

♦  KJ  ♦  (Y  4  (2)  *  (J2-J)  *lk  (3)  •  (J-J1)  )  •  (K-Kl)  J/DJDK 

IP  (K.  GT  .  K  1  .  AND.  J.GT.J1)  HH(K,J)=H4T 

2 10  CCVIIKUE 
2 SO  CONTINUE 

LVAifl-  1 
Jfl= 


DO  300  K=1,KflAX 
DO  280  0=1,JSAX 

IF  (J. IS. 1)  00  TO  26  1 

IF  (LV4B  (K,  J  -  1)  .  GT.O)  LV AE  (K ,J ) =-  1 
2b  1  IP  (XL  (K,  J)  .  E..-999.  .OB.  IL  (K.J)  .  EQ.-999.  )  GO  TO  280 

K  0=  K 

J  fl—  HA  XO  (JB,  J) 

LVAB  (K, J)  =L VABH 
X  (LVABB)  =  XL  |K,  J) 

Y (LVABB) =YL  (K  ,J) 

a (lvabb) =an  (A.j) 

IB  (LVABB) —2 
fl  AT  =  BB  ( K  ,  J) 

IF  (BAT. EQ. 0)  GO  TO  260 

Al2<*-=0.  5*(Xl(K,J-1)  •(XL  (K,  J)  -  YL  (K-  1 ,  J)  )  -XL(K,  J)  *  (XL  (K,J-1) 
a -XL  (K  -  I.  J)  )  ♦XI(X-I.J)  •(IL(K,J-1)-YL(K,J))  ) 

A234-0.  5*  (XL(K.J-I)  *(IL  (K-1.J)  -  XL  (K-1.J-  1)  )  ♦XL  (K-  1,  J)  * 

<t  (YL  ( is  —  1  ,J-1)  -  YL  (K.J-1))  ♦XL(K-  1,  J-  1)  *  (YL  (K  ,  J- 1) -IL  (K-1,J))  ) 

A 2=0.  25*  (XL  (K,J)  +XL  (K.J-1)  *IL  (K-1.J)  »XL  (K-  1,J-  1)  ) 

Y2  =  0.  2d*(YL  (K  ,  J )  *TL  (K,  J-  1)*IL(K-1,J)*XL(K-1,J-1)) 

I L  V  A  a  fl )  =  E  H  C  (MAT) 

I  F  (IP  Co  (BAT)  .  KE.  0)  D  ( LV  AhF.)  =  P.  liO  (BAT)  *  ALPO  (BAT) 
:P<IPOs(flAT).ii£.0)  Ai.  P  (LVABB)  =  ALPO  (BAT) 


A  ii-VALS)  =  A  1  24*A2J4 
Z  (LVABB)  =  D(iVAEK)  •  A  (LVABfl) 

IF  (YC  (SAT) .  BE. 9.)  Y  Y  (  L  V  AF.  K)  =Y  C  ( B  A  T) 
2b0  CONTINUE 
C 

LVABf.  =  LVAEfl*  1 
zeo  CCNTISJE 

300  cost: k j e 

c 

K  fl  A  i  =  A  fl 

J  fl  A  Y  =  J  fl 

•  SITE  (3,  132)  KflAX.JBAX 
132  FOB  BAT  (  IX.  •  r.fl  AX  =  •  .IS,  5X  ,*  Jfl  AY  =  *  .15,/) 


■.--INITIALIZE  (r  L.  IN  CUE  BLOCK 

IP  (I VTYPE.  E_.  0)  GO  TU  450 
LEAD  (  1,  1032)  JD.JU, KU.UZLBu 
»  PI  T  E  (3, 2143)  JH, JO.K 0, 02EBO 

2143  FOB BAT ( li, • J3-* ,12, 3X, • JU=’ ,1 2, 3X , • KU=* ,1 2, 3X. • UZERO=* ,E12. 5) 

C 

AflASS  =0. 

BflA53  =  0. 

B  P  U  =  1  . 
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c 


HKU  1=  1. 

U2INT=UZEBC 


CC  j  ID  J=2, JJ 

IF  (MS  (KU,  J)  .LE.  0  -08.  BB(KU*1,J).  LE.  0 }  GO  TO  310 
ra-bb  |K'i,  j) 

IP  (IPOE  (Hi)  -HE.  0)  BKO=ALPO  (HA) 

BB  =  flfl  (KO*  1,  J) 

IF  (IPOH  (BE)  .ME.  0)  BKU  1  =  ALPG  (S  B) 

ABASS=AB  ASS  *6X0*8 HO  (HA)  *  (XL  (K  U,  J)  -XI  (KO-  1  ,  J;  ) 

BBASS=BBASS*8KU  1*BHO(BB) * (XL(KU*1 ,J) -XL (KO, J)  ) 

3  10  CONTINUE 

IF  (ABASS*BBASS.GI.O..  AN0.IVTYPE.EU.~1)  UZ INT* UZEBQ* BB ASS/ 

J ( Afl  AS  S* EB  AS  S) 

IP (ABA5S*BBASS. GT.O.  . AND.  IVTYPE.EQ.  1)  U 21 NT=U ZE BO* ABASS/ 
t  (ABASS* BBAS5) 

0 

00  325  K=  1 , KB  AX 
DO  325  J-  1, JBAX 
IP (LVAS  (K,J) . LE.O  )  GO  TO  3  20 
Lfl=LVAB  (K,J) 

IP  (K.  GT.  KU.  AMD. IVTIPE. EU. -  1)  I0(LB)=02EBO 

IF(K.IT.KU.  AMD.  IVTYPE.EU.  1)  ID(LH)=UZEHO 

IP  (J.  LE.  JU.  AKD.K.  EU-  KO.  AMD.  J.  GE.  JB)  XD  (LB)  =  OZIKT 

S. 

i  *  V  CONTINUE 
325  CCMT1MUE 

C 

IP  (KCHEK. ME.O)  GO  TO  450 
KCH  EK  =  K  B A  X 

IF  [IVTYPE.EC.  1)  KCH  EK  =  K  0 ♦ 3 

C 

*53  CCMTINUE 

MrIT£(3,125C) 

ZEBO=0. 

DO  470  K  -  1,  K  fl  A  X 
DC  iuO  J=  1,  JBAX 
LB^LVAo  (K  ,  J) 

IP  ( i-.B.  LE.  0)  GO  TO  46  0 
IP (BB (K, J) . GT.O)  GO  To  455 

MBIT  i  (3, 1280)  J,K,Bfl  (K,  J)  ,  LR,  X  (LB)  ,  Y  (LB)  ,  2  E  BO ,  Z  EB  0,  ZEBO  ,  Z  EB  0 

♦  , XO  (LB)  ,  YE  (IK)  , ZERO 
UC  TO  460 

45.  nAT  =  “.n  (K,  J) 
m  =  3. 

IF  (YC  (BAT) .  ME.O)  YYY-YY(LB) 

•  FITE  (3,  1  2  o  0)  J, K.BAI  ,LB, X  (LB) , Y (LB)  , A (LB) , D(LR)  , Z  (LB)  , YY Y, XD  (LB) 
♦ , YD (LB) , E (LB) 

460  CCMTINUE 
470  CONTINUE 
L 

10  10  FCin.-IAT  (£15) 

110.  FChBA7(3F10.3) 

1C  1 1  POBBA1 (UE 10. 3,1 5) 

10  12  Fun. BAT  (4E  10.3,15) 

1  0  3  3  FCfcflAI  (21 5, 4E  10. J,I 5) 

10  32  F03BAT (315, E  10. 3) 

1/  51  FOaBAT  (  IH  1 ,  4i,,J,,4X,,i.,,4X,*H•,,  LV AH • ,  1  IX , • X •  ,  1  IX  , 

♦  IH t ,  1  IX,  In  A,  I  IX,  IriD,  I IX, Ih 2, 7X , 5 HY IELD,  I  OX , 2 HX D ,  10 X , 2H YD ,  I  1 X , 

♦  1  d  E  ) 

1.0  5  PC:(", AT  (415,  5F  I. .0,4  (  1  X ,  E  1  0 .  6)  ) 

C 

a  E  T  U  ?  N 

E  NO 


jUSBuUTINE  3MEEP  (JSIZE, JXI.KXi, XL,YL. BB,LVA8) 
I RPLICIT  BEAL*o  (A-U.O-2) 
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*-  COaaON/T/X (3000)  .If  (3000  J  ,  XD  (3  000}  ,10(3030}  , a  (3000 )  ,  A  (3000) 

♦  , 2(3000] , D (3000) .SIX (3000) ,SYY  (3000) ,S22  (3000)  ,TXY (3000) 

♦  ,TXX  ( 3000)  ,  TYY  (3000)  ,TZZ  (3000) , P{3000) , E (3000)  ,YY  (3000)  , 

♦EW  (3000)  ,EXX  (3000)  ,EYY  (3000)  ,  E2Z  (3000)  ,EXY  (3000)  .ALP  (3000)  , 

♦  IH  (3000) 

COaaON/£vS/EvSrc  (6)  ,  EWSTD  (oj  ,  EtfSTS  £o)  ,  RdO  (6)  ,  1C  (6)  ,  on f’  X 6) 

♦  .CLIN  (o)  ,  C0S-J  (o)  ,T£IQ  (6)  ,SP  (6)  ,  E'OSTG  (6) 
COflBON/GEH/UZEt.O.DT.DTN.DTK.I  JBUHD,  JBAX.JBIN,  KHAX,KHIN,KCHEK 
COBBCN/IND/  NEB  (10) 

coaacN/FSB/  tsr  ( ioj 

COHSCN/POE/  ALP  (3000)  ,SXXS  (3000) ,SYYS  (3000)  ,SZZS(3000) ,TXXS  (3000) 

♦ , TX YS  (3000)  . TZZ S ( 3000) ,TXYS(3  000)  ,PS  ( 3000 ) , EV VS  (3000)  ,EXXS(30  00) 

♦  .  EYYS  (3000)  .EZZS  (3000)  .EX  YS  (3  00  0)  ,EJ  (3000)  ,  ALPO  (6)  ,IPOE  (6)  ,  CTCS  (6) 
O IBEN SION  Xl(KXX.JXX) .YL(KXX, JXX)  ,BB(KXX.JXX)  ,  L  VA  R  ( KX  X,  JX  X) 
DIBFNSICN  XT  EBP  ( 100  )  ,  Y?  EH?  (  10  0)  ,XDT£BP(  100)  ,YDTEBP(100) 

<■  DC  123  1=1,100 

XTEBP  (I)  =0. 

ITEBP  (I )  =0. 

XDTEHF  (I) =0. 

133  X  DTEBF (I) =0. 

e 

C 1  S  «  fl  -  I . 
c 

Du  950  K=  l,  Kfl AX 
DC  920  J=  1.JHAX 

L 

LVAEfi  =  LVAH (K.J) 

IF  (LVABB.LE.0)  GO  TO  780 
C 

V  »  -  0  . 

TXX«i=3. 

T  Y  Y  N  =  0 . 

TZZW  =  0- 
TXI N=  0. 

SIII.=  0. 

SYY  «~0. 

S  ZZ  W  =  0. 

EN=  0. 
py=0. 

0  • 

SPSw=0. 

L 

C - MOHENIUB 

C 

XDKH=XD  (LVAEH) 

Y DNH  =  YD  (LVAES) 

FX=0. 

P  Y=0. 

X  HOfl  =  0. 

A  fl ASS -  3. 

13=LVASB 

u 

C — FIND  THE  COORD.  OP  CELLS  ABOUND  POINT  (K.J) 

C 

CC  360  I*  1,  <1 

u 

C  NBITE  (3,  145  6)  K.J, I 

C  1  4  56  FCRBAT  (  II,  * - '.315) 

DBASS=0 . 

GOTO  (230, 240, 250, 260) , I 

u 

C--I=1,  UPPER  RIGHT  HAND  yUADRANT 
C 

230  IF  (K.Ey.KBAX.OS. J.EQ. JHAX)  GO  TO  360 
IP(BB(K*l,J*1).LE.O)  GO  TO  360 

c 

L  1  =  LV  AR  (£♦ 1,J*1) 

L  2  =  L V  AB  (K , J  *  1 ) 

L4=L VA8  (K ♦ 1 , J) 
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LB  =  L  t 

BAT«BB(K*1,J*1) 

GO  TO  270 

C 

C  —  I=2,UPPEB  LEFT 

C 

240  IF (K. EQ. 1.0B. J. EQ. JBAZ)  GO  TO  360 
IF (BB  (K  ,  J ♦ 1 )  . LE. 0 )  GO  TO  360 

t_ 

L  1=LVAB  (K-1  ,J*  1) 

L2  =  L V  Afi  (K-1 ,  J) 

L4=L?AB  (K  ,  J  ♦  1) 

LB-14 

BAT=BB  (K,  J*  1) 

GO  TO  270 

*_ 

0-1=3,  LOU  EB  LEFT 

250  IF  (K. EQ.  1.0B. J. EU. 1)  GO  TO  360 
IF(flfl (K,J).IE.O)  GO  TO  360 

c 

L1=LVA2  (K—  1  , J- 1 ) 

L2=LVIB  (K, J- 1 ) 

L4=LVAB  (K- 1 ,J) 

IB=L3 

BAT'HS (K, J) 

GO  TO  270 


C - 1=4,  LOUEB  BIGHT 

L 

200  IF (K. EQ. KBAI. OB. J.EQ.  1)  GO  TO  360 
IF(Bfl(K*1,J).LE.0)  GO  TO  360 

L 

L 1=LVAB (K* 1  ,J-1J 
L2  =  LVAB  (K*1,J) 

L4=LVAB  (K,J-  1) 

LB=L2 

fl  AT  =  fl  3 ( K ♦  1, J) 

GO  TO  270 


*70  CONTINUE 

3  00  X02=0. 5*  (X  (LI) ♦  X ( L2) ♦ X (L3) ♦ X (L4)  ) 

Y02  =  0.  S» (Y (II) ♦ Y (L2) *I(L3) *Y (L4) ) 

c 

3  05  A  0=  (X  02 -X  (L 3}  )  *  (Y  U2)  -Y  (L4)  )  *X(L2)  *  (Y  (L3)  *Y  (L4)  -Y02)  *1  (L4) 

♦  • (Y02-Y  (L2) -Y  (L 3) ) 

43=  X  (L4)  *  (Y  (L2)  -Y  (L3)  )-X  (L3)*  (Y  (L2)-Y  (L4)  )  ♦  X  (L2)  *  (Y  (L3)  -I  (L4)  ) 

A  X Y  =  (40* 43) /8. 

AXX=  (Y(L2)-Y(L4))/2. 

4  YY  =  (X(L4)-X(L2)}/2. 

TZZ AX  Y=  0. 

C 

IF (DBASS. NE.O.)  GO  TO  330 
D3ASS=D  (LS)  * 4 X Y 
GO  TO  330 
C 

3  jO  QXX=0. 

£  Y  Y  =  0 . 

2  X  Y  =  0 . 

i. 

c--ST:-4IUS  ABE  POSITIVE  IB  TENSION 
C 

IF  (TBIg  (HAT).Ew-O.  .  OB.  43 .  LE.  0.  I  *  AX  I)  GO  TO  340 

t 

EDXX  =  (  (XD  (L2)  -XD  (L3)  )  •  ( Y  (L2)  -  Y  (L4)  )  -  (ID  (L2)  -XD  (L4))  *  (Y(L2)- 

♦  Y(L3)))/13 

2DYY=-(  (ID (L2)  -  YD  (L3)  )*  (X  (L2)  -X  (L4)  )  -  (ID  (L2) -YD  (L4)  )  *  (X  (L2)  - 

♦  X  ( L3) ) ) /A  3 

EDX I  =  (-  (X D  ( 12)  -XD  (L3)  )  *  (X  (L2)  -X  (L4)  )  ♦  (XD  (L2) -XD  (L 4)  )  *  (X  (L2)  - 

1  I (L3))  ♦ (YD (L2)  -ID(L3)  )  * (Y  <L2) -Y (L4) ) - (YD (L2) -YD (L4) ) *(Y(L2) - 

2  Y  (  L  3)  )  )  /  A  3 
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non  Ci  n  non 


-To  I ABGL  E  Q  STE  ESSES  AE  E  POSITIVE  IS  TESS10S 

CC£F=DSs)BT  (13)  *SP(H1T)  *D  (LH)  *TBIQ  (HAT) 

UXX=CCEF* (2.* EDXX-EDYI) 

2Tr  =  COEF»  (2.* EDYY-EDXX) 

QXY=3.*CCBF*EDXY 

340  FX*EX*(TXX(lfl)*QXX)  *1XX*  (TXY{£8)  *iJXY)  *AYY 

FY-  PI  ♦  (TXY (IH)  *QXY)  *1U  ♦  (TYY  (LH)  *QY  Y)  *111-1221  IT 

1S1SS*1  HASS+DHASS 

KBIT  £  (3  ,  197  6)  l  1,  L2,  E  3,  14,  fl  AT  ,  10,  13  ,  All  ,  1  YY 
U  BIT  E  <i,  179  8)  Q2X. Qtt. Qlt ,*X. 72 ,k*kSS 
C  )  978  POEHAT (515, 45  12. 4) 

C1798  FJ8HA7 (6E  12.4) 

V, 

jc.0  COMISJE 
C 

L  —  tU  >,  AM  V  £  I  . 

L  B- L  V  A3  (F,J) 

IJbABS=IABS  (IJSUKO) 

If  (J.  EJ.  JHAi.  AND.  (1JBABS.  IQ.  l.OE.  IJEA8S.LJ.  5)  )  GU  TO  347 
IF  ( J .  I'*  •  1.1ND.  IJBUNO.SE. -3)  GO  TO  347 
J45  YDh'H=Y0  (LH)  ♦D7N*PY/ABASS 
j  4  7  /HW  =  Y  (LH)  ♦YE1.H*CT 

IF  (K.  £.  .  1  .  A  A I .  ;IJBABS.Ew.4.0E.IJBAr)3.E..5.0h. 

•  IJoABS .  o)  )  GO  TO  357 
I  r  (ft.  £„..K(!AX.  AM.  IJBABS.  Zj.  4)  ,G  To  357 

c 

jSp  X  JNn  =  XJ  (IB)  ♦0Ti.*PX/AH13S 

1 

3  c  7  X  Si.  =.(  (LH)  ♦  X-  Sit*DI 

IF  (K3  (ft, J)- Ew.O)  GC  TO  750 

L  -  -  N  E  •  kk  r.  k  AM  VC  I.  FCE  CELL  ft  ,  J 

A  I2*4=XTEHP  (J-  tj  •  (YMi-YTEBP  { J)  )  -XM*  (YTEHP  (J-  I)  -  Y  T  E  B  ?  ( J )  )  * 

1  ITIHf  (J)  •  (YTEHP  (J-  1)  -YHS) 

L  H  -  L  V  AE  IK, J) 

LRK-LVAB  (A-  1,J-  1) 

LAB  =  LVAB(K#J-1) 

LHJMVAo  (A-  I.  J) 

A  -  3-*=  XT  EH?  (J-  1)  *  (YTEHP  (J)  -YKH  JH)  ♦1TEHP  ( J)  »  (YaH  >H-  YTEHP  (J-  1)  ) 
I  XKHJK*  (YTIBB(J-I)  -YTEHP  (J)) 

c 

423  lt=3.  5*  (A  Il4»A234; 

IF  (Ab.3T.J-  i  Jo  To  425 

■  SITE  (3  /  1013)  V  ,  J 

10  10  FOfi.H  AT  (•'•*.  *»  k.  ■■  rfW  «W  WBUHarr.r.  f  f  BRf..-  £2  5  1  ,  21  5) 

4  I ‘  C-hTISJE 

I‘  =  7  (I  “■>/!■ 

■  RITE|j,I5l7j  r.#J#AI24,A234,Alr,Dn 

2  15  6'  F  L  f  r 1 T  (  1  X , 2  I  -  ,  - f  15.  4 ) 


C- -CM  F  ETE  ST  HAIM 
o 

II  A  =  31/  (  A  «  *  A  (Lflj  ) 

X  1  J  -  (I  (IF)  •XNW)/2-  -X8.1H 

X  H  -  (XTEFP  (J  -  1)  *  X  (LKfl)  - 1 7  E  fl  ?  (  J  <  “X  i  L  fl  J  )  )  /  2  . 

T  ’  1  -  -  (  Y  ( I  H)  ♦  Y  N  p  )  /  2 .  -  Y  HH  R 

Y  :•-£=  (YTEHP  jJ  -  1)  ♦  1  (LrvHj  -  Y  T  E  HP  ( J  j  -  Y  (  LH  J  j  )  /  c 

1  I.;  1  J  =  »  J'.i  -  I.  ( L  H  A  ) 

XIMMMTEHt  ( J  -  1 J  -  X  L  T  EH  P  ( J  ) 
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n  r 


YUH13=YDkH-YD  (Lflk) 

YDH <4 2  =  YDTE.tr  (J-  1)  -YDTEMP(J) 

— DEPIHt  ccoed.  cf  ceil 

X  1=0.5*  (X  (L  6)  *iN») 

I  2=0.  5*  (XT EBP  (J)  «X  (LflJ)  ) 

X 3=  XHBfl 

I  <4=0.  5*  (XTEfl?  (J-  I)  *X  (LKfl)  ) 


Y  1=0.  5»  (Y  (LF)  ♦Y  MH ) 

Y  2=0.  5* ( Y  T  E  B  P  (J)  ♦Y (LB  J)  ) 

Y  3=  YFSB 

Y4  =  0.  5»  (YTEflP  (J-  1]  ♦  Y  (LKfl)  ) 

C 

C--S„i)AEED  CF  VECT0B3 

C 

XHAG4  3=  (X<4-X3J**2*(Y<4-Y3)*»2 
v  HA  G4  1=  (X4-X  1)  **2*  IY4-Y  1)  **2 
X  HA  G  1  2=  (I1-X2) *  *2  ♦  ( Y l-Y  2) *  *  2 
XflAo2  3=  (X2-X3) **2*  (Y2-I3) **2 

c 

C--DOI  FRCD 

c 

D4J2=-  (  (X4-X3) *  (I3-X2)  *  (Y4-Y3)  *  (13-12)  ) 

DJ2 l  =  - (  (X3-X2) *  (I2-X1) ♦  (Y3-Y2) •  (I2-Y  1)  ) 

Dl4  3=-((Xt-X4)*  (X4-X3)  ♦  (Y  1-14)  *  (14-  T  3)  ) 

0  2  I  4  =  -  (  (X2-X  1)  *  (X  I- 14)  ♦  (Y2-I1 )  *  (Y  1-14)  ) 

4. 

C--  CLECr.  TC  SEE  IF  PHOJECTIOB  LIES  IHSIDE  CELL 

IF  i  D432. LE.O.)  D432=  0. 

I  F  l  D32I.LE.O.)  D32  1=0. 

IF)  D214.LE.0.)  D2  1  4=  0. 

IF(  D143.LE.0.)  D  14  3=0. 

c 

D<*  32=  C4  32«»  2 
Di2  1  =  D3 21**2 
D2  14  =  22  I4»*2 
D  143=2  I43»*2 

C 

<.--FIHD  TEE  fllMIHUfl  DISTANCE 

DELX  =  D.tIN  1  (IB  AG  43-D43 2/XH AG 23  , 

1  XH AG  23-D432/XH AG43 , 

2  XH  AG  2 3-D 32  1/XK AG  12  , 

3  XflAG  12— D32  1/XHAG23  , 

u  XflAG  12-D2 1 4/XHAG41 , 

a  XHAG4  1-D2 1 4/XflAG 12  , 

c  XHAG4 1-D 14  3/XHAG43 , 

7  XflAG 43-D  14 3/ XflAG 41  ) 

£VUL=2.  •  (D  (Lfl)  -DU)  /  (D  (Lfl)  ♦  OH) 

EIIil=  DT  A  •  (X0H42*IH13-IH42*IDH  13) 

EYYH=-DTA*(Y0H42*XH  1  3-XH4  2  *ID  H  1  3) 

EXT  H  =  0. 5*DTA* ( I DH42 *T H I 3-TH42 *YDH 13-XDH 42*XK 1 3*XH42*XDH 13) 
EZZH=EVCL-EXXa-EIia 

A L F  A  =  0 . 5  *DT  A  *  (-IDH4  2»ia 13* I H4 2* YD H 1 3-XDH4 2 * XH 1 3 +1 H4 2* XDH 1 3) 
430  flAT=flfl(K,J) 

X- 

c  kEITE <3 ,3562)  J, K, D ELX , ETOL , E XXH , EY IH, EZZ a, EX YH , A LF A 

CJ‘f.  FCfiflAT(2I5,7E12.4) 


C -- VI SCGSI T Y 
SFF.=  1. 

If  (IPOS  (HATj.NE.O)  SPa=  1.  ♦  (  1.  -  ALP  (Lfl)  )  /  ( 1.- ALPO  (HAT)  )  *  ( 
♦  CTCS  (BAT) - I-) 

L 

DELD=Dk-D  (LF; 

IF (DELD. GT. C. ) 
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n  n 


♦  J^DELD/DT*  (SPB*SP  (BAT)  •CLIN  (BAT)  *DSQfcT  ( AH)  »CySQ  (BAT) 
♦•AN*D£LD/DN/DT) 

C  IF  (DELD.LT.  0.  AND.  IPOB  (BAT).  NE.O) 

C  *Q=D£LD/DT» (SPR*SP (BAT) *CL I N  (H  AT) *DSQRT  (AH) ♦  CQSQ  (BAT) 

C  ♦*AW*DELD/DH/DT) 

C 

C--E3TIHATE  INTERNAL  ENERGY 
C 

D  ELZ  =  (SIX  (LB)  •  EXX H* S Y  Y  (LB)  *EYYH*SZZ  (LB)  *EZZH*2.  *TXY  (LB)  *EXYH)  /Z  N 
£«- E  (LB)  *DE1Z  -  (P(Lfl)  ♦.)  •  (  I.  /DH-  1.  /D  (LB)  ) 

EAV J=EVOL/3. 

5  £T  A=  2 .  •  T X Y  (LB)  ‘ALFA 
C 

IF  (IFCS  (BAT) . EJ.O)  GO  TO  60  1 

CALL  FGREL (EH, E VOL, EXXH,EITH, EZZ H , E X Y H . S X X H , S Y Y W# SZ Z H , T X Y H , PH , 
Eir.E'IIA.ALFA.LB.BAT.K.J) 

1C  6  0  0 
c 
c 

C--ELASTI«-  BA1ERIIL 
C 

60  1  <-  ALL  E  *  SI  (£«,DH,P.,BAT) 

C 

c 

SXXk  -  SIX  (LB)  *2.  *U8U  (BAT)  •  (EXXH-EAVG)  ‘BETA 
S  Y  YN=  SY  Y  (LB)  ♦*.  *0BU  (fl  AT) *  (£YY  H- EAYG) -bETA 
SLZ«=SLL(Lfl)*2.  »UBU  (BAT)  *  (EZZH-EA VG) 

C 

7XY  h  =  TX  Y  (LB)  *UMU  (BAT)  *EX  Yd*  (SYY  (LBJ  -SXX(LB)  )  *ALFA 


C 

IF  (YC  ('.AT).  I  E  •  0  .  )  GO  TO  bOO 

SJ2  =  SIX ■••2  *S YY N  •  *  2  ♦  SZZ  H  *  *  2  ♦  2 .*TXYn»»2 

YYY=0.66o6t7*YY  (LB)  »«2 

I F  (SJ2. LE. Y YY)  Go  TO  600 

CI=)jvol (YY  Y/SJ2) 

SXXH=CY*SXX > 

SYYH*CY*3YYk 

SZZN=CY*SZZN 

TXYH=CY«TXYk 

c 

C— ADJUST  IBTEBNAI  ENERGY 

C 

o  Ou  Ek=£  (LB)  *0.  5*  ((SIX (LB) ♦SXXN)*SXXH* (SYY (LB) ♦ SI Y H) * EY YH ♦ (SZ Z ( LB )  ♦ 
1  SZZN)  *EZZa-»2.*  (TXY  (LB)  ♦TIIW)  *EITH) /DH-  (  ( P ( LB) ♦ PH ) /2 . *Q) *  ( 1 . / 
20M-  1./D (LB)  ) 

o 

--  COBPUTE  TOTAL  STRESS 

620  TXXH=SXXB-PH-3 
TYYH=SIIH-PH_y 
TZZH=SZZH-PH-Q 

C 

C - SEPARATION  OF  IBPACT  PLANE 

C 

IF  (HER  (BAT)  -HE.  5.  AND.  NPH  (BAT)  .  N  E.  6)  GO  TO  690 
IF  (NFE  (BAT)  .EQ.6)  GO  TO  650 
IF  (TXXH. LE. TSB(BAT)  )  GO  TO  69  0 

C 

P  1=TXXH»EwSTC  (BAT)/  (EOSTC  (BAT)  ♦  I.  3333*UHU  (HAT)  ) 

DSL  =  1.3333*TXXN*OBU  (HAT)  /  (EQSTC  (BAT)  ♦  1. 3333*UBU  (HAT)  ) 

C 

T11N  = IXXH-OSX-P 1 
TYYH-TYYH-P  l*DSX/2. 

TZZN=TZZH-P  )♦ DS X/2. 

r 

v. 

FH^-(TXXH*TIYH«TZZH)/3. 

TXYH^C. 

GO  TC  600 
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cSO  If(TXYN.LE.ISB (BIT) )  GO  TO  690 

t 

P 1=TYIW *EQS1C  (BAT) /  (EQSTC (BAT) ♦ U 3333*UBU (BAT) ) 

DSY  =  I.  3 333* T I IV *0 BO  (B AT)  /  {EQSTC  (BAT)  *  t.  33  33*0  BO  (BAT)  ) 

TXX  W=TXXi-P  WDST/a. 

TYYN=TIIN-P  1-DSY 
TZZW=TZZW-P  1*DSY/2- 

c 

PN=-  (IXXN*TYYNMZZW)/3. 

1  *y  «=  0- 

c 

b  60  IT  (18  (LB) -EC-  1)  GO  TO  690 
18 (LB) *  1 

C 

WHITE  (3  ,  1680)  K,J 

1680  PCBHAT  (  IX,  • IEPABATIOW  AT  CELL  K,J  =',2I4) 

C 

C--COSPUTE  SOUND  SPEED  AND  TIB.  *~t' 

c 

o90  EBOD=0. 

SPSQ=SP  (BAT) **2 

If (IPCB  (BAT) . HE. 0)  SPSQ=SPSQ* SPB*SPB 
If  (DABS  (DW-C(LBJ).LT.  1.  B-8)  GO  TO  700 

EHOD=PW/  (DN/BHO  (HAT)  -  1.  )  *2.  »Q*DN/  (Dtl-D  (LB)  )  ♦  I  .  33*0H0  (BAT) 
SPSU=DHAI1  (EBOD/D  (LB)  ,  .3*SPSQ) 

7  00  D1S  2= DELI/S P SC 

If  (DTS2-GE.DT SOB)  GO  TO  750 

KT=K 

J1=J 

DEL XT=DELX 
DTSJT=DTSC 
SPSJT=SPSQ 
DTScfl*3TS0 

c 

7;0  CONTINUE 
7.0  CONTINUE 

V. 

If  (K. EQ.  1)  GO  TC  790 
If  (J. EC.  1)  GO  TC  785 
L8N=L  VAB  (K-  1,  J-  1) 

If(LBH.LE.O)  GO  TO  785 
I  (LBN) «XKBJB 
Y  (LBN)  =  YKHJB 

7  d  5  CONTINUE 

LHJ=L?AB  (K-  1,  J) 

If  (LBJ. LE.O)  GO  TO  790 
XHBH=  (X  (LBJ) ♦XTEBP(J) )/2. 

YUBB=  (Y  (LBJ)  *  YTEBP  ( J)  )/2- 
XKBJB-XTEBP  (J) 

YKBJfl=YTEBP  (J) 

ID  (LBJ)  =XDTEBP ( J) 

YD  (LBJ)  -YDTEBP  (J) 

C 

790  If (LVABB.LE.O)  GO  TO  920 
XTEBF  (J) =XNV 
YTEBP  (J) =YNS 
XDTEBP(J) =X £  N  8 
YOTEBP  (J) *YDNH 
If  (Bfl  (K„J) . EC.O)  GO  TO  800 
c 

795  LB-LVASM 
0 (LB) -D W 
E  (LB)  *£N 
SIX  (LB)  =  SXX  k 
STY  (LB)  =  S  Y  Y  k 
SZZ  (LB)  -SZZk 
TXY  (LB)  *TXYk 
TXX  (LB)  -T  XX  k 
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TIT  (LB)  =TYYH 
TZZ  (LB)  -TZZ h 
P (LB) =  PH 

BIT  (LB)  =  Ef?  (LB)  *EVOL 
EXi  (LB)  =  EXX  (LB)  ♦  EXXH 
Bit  (LB)  =EII  (LB)  *EIIH 
EZZ  (LB)  *EZZ  (LB)  ♦  EZZH 
EXT  (LB) =EXY  (LB) *EXYH 
ILF  (LB)  =XLF  (LB)  «ALP  A 

800  COHTIHUE 
920  COHTIHUE 

I v  (K. IQ. 1)  GO  TO  940 
L  BJ  =L  TAB (R- 1,JBAX) 

IF(LBJ.LE.O)  GO  TO  940 
X  (LB  J )  =  X  KB  J  B 
Y  (LBJ)=IKBJE 

KK-K 

IF  (K. LT. KCH EK.OB.K. EQ. KBAX)  GO  TO  950 

DO  945  J  =  I, JBAX 
LSJ -L  VAi  (K-  1,  J) 

IF  (LBJ. LE.O)  GO  TO  945 

IF  (DABS  (XD(LBJ)  )  .  GT ..  1E-3.0B.  DABS(YD(LMJ)  )  .  GT.  .  IE-3)  GO  TO  94  8 
945  CCHTI HUE  ' 

GC  TC  960 

948  KCHEK=BIHO(K*  1.KBAX) 

GO  TC  960 
950  CCIITI  HUE 
960  CCHTI HU  E 

DO  960  J- 1, JB AX 
LBJ=LVAB  (KK  .J) 

IF  (LBJ, LE.O)  GO  TO  980 
1  (LBJ) - XTSflt  (J) 

I  (LBJ)  =YTEBF  (J) 

10  (LBJ)  = X  DT E H P ( J) 

YC  (LBJ)  - Y  DT  EBP  (  J) 

CC HI I  HUE 

DT.  =  DSJET  (DIS^B)  •  (  1.-3.  *TEI  J(  ) 
letuen 

E  HD 


SUB BOUT  I  HE  iCST  (£,D,P,fl) 
IB ELICIT  KEAL*8  (A-H,0-Z) 


CCKBOS/  E  .S/l^STC  (6)  ,  E  y)STD  (6)  ,  EjSTS  (b)  ,EH0  (6) 
*  ,  CLI  b  (6)  ,CUSJ  (t>)  ,TiIJ  (6)  ,SP  (6)  ,  EJSTG  (6) 
EBU=D/EHu  (B) -  1. 


YC  (6)  ,  UBU  (6) 


PH=EB0»  (E.S1C  (fl  )  ♦  EBU*  (EJSTD  (B)  ♦  EBU*  E'JSTS  (  B)  )  ) 

?=PB*  (1  .-0.  5*E2STG(H)  *(  1.  -  E  HO  (B)/D)  )  ♦EjSTG(fi)  »RH0  (B)  *E 
EEIUBH 


E  HD 


C 


SUBROUTINE  FOBEL(DH.E?OL,EXXH.EYYil.  EZ  Z  H  ,  E  X  Y  H,  SX  XW  ,  S  Y  Y  H, 
♦3ZLH,TXYN,PH,EN,BETA,ALFA,Lfl,flAT,K,J) 

I BPLICIT  BEIL*8  (A-H.O-Z) 

CCBBGH/T/X(3000) ,1(3000)  ,XD  (3000)  , Y D  (3000 ) , fl ( 3000 ), A  (  30  00 ) 

♦ , Z ( 3 0 00 ) ,D(3000),SXX(30  00) ,SYY  (3000) ,SZZ(3000)  ,TXY(3000) 

♦  ,TXX  (30  00) , TIT (3000) , TZZ  (3000 ) , P (  3000) , E (30 00)  , IT (3  000)  . 
♦£W  (3000)  ,  EXX  (3000)  ,  EYY  (3000)  , EZZ  (3000)  ,EXY  (3000)  ,  ALF  (3000) 
»i a (3000) 


CCBBON/EiS/EQSTC  (6)  ,  EQSTD (6)  ,  EJSTS (6)  ,BUO (6) 
♦  ,CLIM  (6)  ,CUSQ  (6)  ,  TBI J  (6)  ,SP  (6)  ,  EQSTG  (6) 


YC (6)  , 0B0  (fa) 


COBBOH/POE/ 

♦  ,  TITS  (3  000)  , 
♦ , EY YS  (3000)  , 


ALP (3000)  ,S  XXS  (30  00) , ST  IS (300  0)  ,SZZS ( 3000) , TX XS (3000) 
TZZ S (3000)  ,TXIS (3000)  , PS  (3000) , EVVS (3000)  , EXX SI 30 00) 
EZZS(3000)  ,EXTS  (3000)  ,  E J  (  3000 )  ,ALPO(6),IPOR  (6)  ,CTCS(6) 
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—  — ALP2,  ALPA  AT  PB  EV 10 US  TIflE  PTO,  PBESSUBE  AT  PREVIOUS  TIBE 
A  L  PP ,  PTP  ,  ABE  VABI  ABLES  FOB  ITEBATIOti  PROCESS 

ALPy= ALP  (LB) 

k  BITE  (3,  123)  K, J, LB. ALP  (LB)  ,P  (LB) 

123  POBHAT  (IX,'  EEPOBE'.aiS.SX.IEU.  5) 

114  F0B8AT (IX, 'AFTER*. 315. 5X.5E  12.5) 

A LPP  =  AL  P *' 

PTO  =  P  (LB) 

PTP-P10 

HH  =  I.  ♦  (  1.  -  ALP Q)  /  (  1.  -ALPO(BAT)  )  •  (CTCS  (BAT)  -1  .) 

11=  1 

10  BHOS=DH/ALPf 

CALL  EQST  (EH.BHCS.PSS.BAT) 

PTT  =  A  LPP  *PS  £ 

IF (DABS  (PTP-PTT).LE.O. 1)  GOTO  100 

- THE  P  HOC  ED  UR  E  TO  FIND  N  EH  PTP 

A  k  U  =  E HOS/ EH C  (BAT) 

DkPP-  (l./HH/HH-  1.J/E0STC  (BAT) 

DAPT=-1.  /  (ECSTC  (BAT) *£0510 (BAT) * (ANO- 1.  ) *  (ABU* I.)  ) 

PTP=  (PTP*DAPP-PII*DAPT) / <D APP-DAPT) 

ALPP=  ALPC*  (FTP-PTO)  /ECSTC  (BAT)  •  (U/HU/HH-  1.  ) 
k  HITE  (3, 100C)  II.LB.BHOS, ALPd.ALPP.PTT.PTP.HH 
1000  FOBBAT ( 21 5, EE  12.5) 

11=11*1 

IF (II. GT. 20)  CALL  EXIT 

GC  TC  10 
100  CONTINUE 

- NEGATIVE  PRESURE  AT  UNDISTURBED  CELL  IS  FOBCED  TO  GO 

LACK  TO  OEIGIkAL  STATE 

IF  (ALFP. LT.  ALPO  (BAT))  ALPP= AL PO  (B AT) 

0 

C - N  E *  ALPA  AND  FBESSUBE  V£i E  DECIDED.  NOW  GOING  FOB  DEVIATOBIC 

C  STRESS 
C 

DALP=ALPP-AIPJ 
E A VG=  EV  CL/ 3. 

ALVT=  1.  *DALE/EVOL/ALPP 
IF  (ALVT.LE. 0. )  GO  TO  2*5 
IF  (ALVT.GE.  I.  )  ALVT=  1. 

GC  TC  247 
^49  AL  VT=  0. 

247  ETXX=EXXH-EAVG 
ETY Y=EY YH-EIVG 
ETZZ=EZ  ZH-EAVG 
C 

ESXX= ALVT»E1IX 
Z5YY=ALVT*E1YY 
ESZZ= ALVT»E1ZZ 
ESX 1=  AL  VT*E  XY  H 

0SXX=2.  *UBU  (SAT) *ESXX*BETA 
GGY  Y=2. *UBU  (BAT) *ESYY-BE7A 
0  3ZZ=  2.  *UBU  (BAT) *ESZZ 

D3XY  =  2.  *  U  BU  (BAT)  *ESXY*(SYIS(LB)  -SXXS  (LB) )  *ALFA 

SXXS  (LB) =SXXS  (LB)  *DSXY 
SY  YS (LB) =SY  YS  (LB)  *DSYY 
SZZS  (LB) =SZZS  (LB) *DSZZ 
7XY  S (LB) =  TA  IS  (LB)  *DSYY 

GXA  k=  ALPP*S  XIS  (LB) *1.0 
SYY  H  =  ALPP*SYYS(Lfl)  •  I. C 
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szzn=alpp»szzs(LH)  •  i.  o 

TXY  N=  ALPP*T  XYS  (LB) • t. 0 
C  RR1TE  (3,  1457)  SXXU,SYYi4,SZZW,ALVT 

CH57  FOR  BAT  (  IX,  *  CE  V  •  ,  4  E  1  5.  5) 

C 

c  :k  cut 

C  AL !-  ? 

C  ALPC(WAT) 

c  ej (La) 

C  ALP. 

CL. 

C  FTP 
C  R HO  (MAT) 

C 

CALL  EL IPT ( ALPP , ALPG (HAT) ,  EJ ( LH ) , ALPQ,DV,PTP,P(LH), BE  0  (HA  T) , 
♦ECCTC  (WAT)  .EQSTD(HAT)  ,S XXV , SY YV , SZZV , TX YV , K , J) 

C 

F  V=  PT  P 

ALP  (LW)  =  ALPQ 
PS  (LH)  =PV/ALPQ 
SXXS  (LB) =SXXV/ALPW 
SYYS  (LH)  =SYYV/ALPQ 
SZZS  (LH ) =SZZV/ALPQ 
TI Y S  (LB) =TXYV/ALPQ 
C 

243  BETUB!. 

END 

C 

C 

SUBROUTINE  ELIPT (ALFO, ALPP, £J P, ALFA  1 , RHO, ? P.PBP,RHOS, IP , B P, SXXV 
♦  , S Y  Y V ,S  ZZM ,TXYV,K,J) 

I HPLICIT  BEAL*8  (A-H.O-Z) 

C 

C - £J , AL  FA  1  AT  PREVIOUS  STEP  FOE  STARTING  POINT 

C 

AAA-  1.00 

y=2.  o:06 

BET  A-  6. 

C - PUT  PREVIOUS  ALFA  FOR  CHECKING  PLASTIC  RANGE 

ALPA=ALFA 1 
EJ=  EJ  p 
P=PP 

QC=  1. 000 
A  H=  2. D00 
c=. ocoo 
CP  =  . 75 
CK=  1.  00  006 

C - CURRENT  DEVIATORIC  TERH 

SUH=SXXV*SXXV*SYYV*STYV*SZZV»SZZV*2.*TXYV*TXYV 
EJO=DSQRT ( 1 . 5  *S  U  H) 

E J 1  =  E  JO 
11=  1 
KK=  1 

<_ 

C - CAROLL-H0LT' S  PEESSUBE  EQUATION 

C 

R  1  =  2. COO/3. COO 

c 

go  tc  no 

■JO  11=2 

CCC  IF (PP.GT.PRF  )  GO  TO  100 

c 

C - PLASTIC  UNLOADING 

C 

CCC  CALL  PL  AST  (P  1 , Y  1,QO,AH,C,OK,PBP,PP, EJP, EJO, EJ 1) 

CeC  KK= 2 

CCC  GC  TC  400 

c 

i. - STARTED  WITH  FBEVLOUS  ALFA 

C 

100  b  2=  R 1 *B  El  A 

R J=AAA- ALFA 


RESULT  FEOW  ELASTIC 
INITIAL  POROSITY 

3J  2*  *. 5  AT  PREVIOUS  TIBE  FOR  STARTING  POINT 
PREVIOUS  ALPA  FOR  STARTING  POINT 
CUBfiENT  TOTAL  DENSITY 
:  RESULT  FEOW  ELASTir 
INITIAL  SOLID  DENSITY 
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P 1=  Y/BETA* (  1. D00/B3**B2-CP/ (AAA-ALPP)  **B2) 

DP1  = I.*B1*T/S3»*(B2*1.D00) 

DDP  1=  1-  *B  1*Y*  (8  2*  1.D00)  /B  3*  *  (82*2.  DOO) 
c 

C - PBESSUBE-ii.fi  B  ELiTIOK 

c 

B  4=  8 HO/ ALFA /BUG  S 
B  5=  B4-  1 . COO 

c 

if  (II. £(..  1)  GO  TO  799 
P=ALFi*  (AP* 5S*BP*Bi*B5) 

DP=-AP-8P*  (B4*R4-  1.  ) 

DDP=2 ,D00*bF/ALFA*fc4*B4 
c 

C - DEVIATOBIC  TEEH  AND  ALFA 

t 

799  Y  l=yG*fi 3**A H*P 1 

DY  J  =  wC»  (-ifl  *R  3*  *  (AH  -  1. DOO)  *  PUB  3** AH*  DP  1) 

CDY  1=  JO*  (AH*  (AH  -  1.  DOO)  *B3**  (AH-2.  DOO)  *P  1-2.  000*  AH*R3**  (AH- l .D 
♦  00) *  DP  I  ♦B3**Afl*DDP1) 


- A.  AND  b 

E6=  1.  O00-C 
57=  1.  E00*C 

i=2. D00*P-B6*P UOK 
DA=2.D00*DP-Bo*DP  1 
DDA  =  2.D30*DEP-f!6*DDP  1 

E?  =B  7  •  P  1  *OK 

IF  (II. Ej.  I)  GCOH=A*A/B/B*EJO*£JC/Y  1/Y  1-  1. 

IF  (II. EC.  1.  ABC.  CCOH.GT.  0.  )  GO  TO  90 
0  WHITE  (3.  120  C)  GCOH, A, B, EJO, Y 1 

Cl. 00  FOBHIT  (  U,  •  GCOH  •  ,5E  15.  5) 

IF  (II. El..  1.  ABC.  GCOH. EE.  0.  )  GO  TO  30  0 

V 

CS=B7*D?1 

DEB=67*DDP1 

V. 

C - YIELD  FUNCTION 

C 

G=A*A/B/B*EJ*EJ/Y  1/Y 1- 1. D00 

Gi=2. DOO*  (A*DA-A*A* DB/B) / (B*B) -EJ*EJ*DY  1/Y 1/Y  l/Y 1 
G  A A- D  A*  DA  *3 • DOO  • A *A • DB*DB/B/B - 4 . DOO *A*D A* DB/B 
G AA  =  (GAA*A*EDA-A*A*DDB/B) /B/B 

GAA=2.D30*(GAA*  EJ*EJ/YI/Y 1/ Y 1  * ( 3. D00*DY  1* DY  1/ Y 1 -DDY 1) ) 

c 

G  J-  2  .  D00*EJ/Y  1/Y  1 

GJJ  =  2 . DOO/Y  1/Y 1 

G AJ  =  — 4.  DOO*  E J/Y  1  *DY  1/Y  1/Y  1 

C 

C -  ELEHENT  CF  JACOBIAN 

r 

AJ  I  1=2.  DOO*  <GJ ♦  (ALFA-ALFO)  *GAJ-  (EJ-EJO)  *GAA) 

A  J  1  /=  2.  DOO*  (  (ALFA-ALFO)  *GJJ  GA-  (EJ-EJO)  *GAJ) 

A J2  1  =  G A 
A J2  2  =  GJ 

C 

D  ET  =  A J  1  l*AJ22-AJ12*AJ21 

C 

C - ANOTHER  NONLINEAE  EO.  WITH  YIELD  FUNCTION 

C 

F= 2. DOO •  (ALFA-ALFO)  *GJ- 2. 000*  (EJ-EJO)  *G A 
C 

C - NEW  ALFA  AND  DEV.  TEBH 

C 

ALFAA=AEFA-  ( A J2 2* F- AJ  1 2 *G) /DE T 
E J  1  =  E J-  (-AJ  2  1*F*AJ  11  *G) /DET 
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nor 


- EJJ  FhO  H  Y I  EL  I  FUNCTION  DIRECTLY 

EJJ=2.D00*DSJRT  (DABS  (  (P  1-P)  •  (P*OK*C*P  1)  )  ) 

EJJ=EJJ/(f  1*CK*C*P1)  * Y  1 

IF  (II. GE. 45)  WRITE  (3,  10 00 ) K , J , 1 1 , ALFA ,R HO . E J, EJ J, P, P 1 , F , G 
1  3 DO  FCFHAK3I5,  8E12.5) 

C 

EE3  =  1.0  C-08 

STOC=  DA  fcS  ((EJ  1-  EJ)  /EJ  1) 

STCD=DAES  (  (ALFA A- AL  FA )  /  ALF  A  A) 

IF  (5TCC.LE.  ERE.  AND. STOD.LE.  ERR)  GO  TO  400 

£J=  E J  1 
ALFA=ALFAA 

FF=F 

GG  O 

II-II*1 

IF  (II. GE. 50)  WRITE  (3,  10  0  1) II. RHO,  STCC , STOD 
IF  (II. GE. 50)  CALL  EXIT 

1001  FORMAT (5X, • NEWTCN  METHOD  DOES  NOT  C  ON  YE  EG  E  •  ,  15 , 3E  15 . 5 ) 

GC  TC  100 

400  CONTINUE 

If (KK. RE. 2)  PP=P 
EEJ  J  =  EJ  1/EJC 
SXXU=SXXN*EEJJ 
SYYW=SYYW*EEJJ 
SZZW=SZ2W*ZEJJ 
TXYW=TXYW*2E JJ 

300  IF  (II. Et,.  1)  ALF A  1=ALFO 
IF (II. NE.  1)  ALFA 1=A  LFA 
Cc  IF (ALFA 1.LE.ALPP)  ALFA 1=ALPP 

EJP=EJ 1 

s z:ucb 


SUBROUT INE  FL AST  ( P  1 , Y  1,  JO,  A M,  C , CK , PEP . ? , E J, EJ O , EJ J) 


IMPLICIT  EEAL*d  (A-H.O-Z) 

PK=OK»C*P 1 
FP=. 5* (P  1  ♦ P  h ) 

P M=  .  5*  (  F  1  - P K ) 

- P.-.F.EJ  ARE  PREVIOUS  POINT  FOT  STABTING  POINTS 

11=  1 
PU  1=PRP 
EJ  1  =  EJ 

- JACOBIAN 

10  AJ  I  1=2.  ♦  (PU  1-PM) /PP/PP 
AJ  12=  2.  *  E J 1/Y  1/Y  1 

AJ2  1  =  4.  *  (EJ  1/Y  1/Y  1-  (EJ 1-EJOj /PP/PP) 
aj::  =  4.  *  ( (?ui-p)  /  y  i/Y  i-  (?u  i-pti)  /pp/p?) 

D  ET=  A  J 1  1*AJ22-AJ12*AJ21 

G'G=  (FU1-PM)  »  (PU  1-PM) /  PP/PP*  EJ  1*EJ  1/Y1/Y  1-  1. 

FF=4. *  (EJ  1*  (PU1-PJ/Y 1/Y  I- (PU1-PH)  *  (EJ 1-EJO) /PP/PP) 

P02  =  -  ( A J  22*  GG- A J 12*  FF ) /DET*  PU  1 
EJ2=  E J 1  - (AJ  11*FF-AJ2 1 *GG)  /DB T 


EAB=  I , D *0 2 

STOC  =  DA  BS  (P02-PU 1) 

STCD=CAES  (EJ2-EJ 1) 


•  / 


\  f  .  l k.  i ii i* .  AND.  j’Oi> . i.i; .  Kii ;  ;.i  tu  jo 

* .•  i t 1 1 ,  i  50  i)  t»yj ,  rv i  ( r  i ,  v(.,  i<;>  s>.i 

i  ;  )..■ 

£.1  I  -  t.J  J 
: ;  - 1 :  •  i 

W  ft  I T  K  (  i,  11  II) 

ii  r>h ~.n  i  u  ,  •  Mirftth  a:.:;c  tasvsiiji. 

I.  i  *  I  .  #  _  .  j  U )  ciLL  i.  XII 
I  y  ( 1 1 .  LT.  :ju  J  .io  TU  I  i 
.  ■>  i’  -  ;>  u  j 
L J J  =  £J  J 

,i.k. :  m  i. 
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